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Mountain pine beetle (Dendroctonus ponderosae; MPB) population has existed at 
endemic levels in the pine forests of western North America for centuries, but in recent 
decades it grew to epidemic levels and outbroke over extensive areas from British 
Columba in Canada to New Mexico in the United States. The current MPB outbreaks 
have impacted large expanses of lodgepole and ponderosa pine forests, reduced their 
ability to act as carbon sinks, altered wildfire hazards, affected wildlife populations, 
changed regional climate, modified local surface energy balance and water quality. Those 
effects are predicted to increase as a consequence of the direct and indirect effects of 
climate changes. Despite severe impacts of MPB, substantial unknowns and uncertainties 
still exist about its historical and current spatial-temporal patterns, future potential 
distributions, disturbance regime characteristics, ways of interaction with other major 
disturbance events, and impacts on forest resilience mechanisms. 

In this dissertation, I first explored the potential of medium resolution satellite imagery in 
mapping the chronic insect disturbance in the Southern Rocky Mountains Ecoregion. A 
forest-growth trend analysis method that integrates temporal trajectories in Landsat 
images and decision tree techniques was introduced to derive annual forest disturbance 
maps over a period of one decade. This workflow is able to capture the disturbance 
events as represented by spectral-temporal segments after the removal of observational 
noises from temporal trajectories in Landsat images, and efficiently recognizes and 
attributes events based on the characteristics of the segments. Higher overall accuracy 
(OA) was achieved when compared with the traditional single-date classifications, and a 
smaller number of training sample units is required compared with maximum likelihood 
and random forest classifiers. 


To test the feasibility of the trajectory-based approach at broader scales, I advanced this 
method by replacing the decision tree based semi-automatic event labeling procedure 
with an automatic attribution step via random forest, which was run on a set of segment 
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features containing information on spatial-temporal neighborhoods. Meanwhile, I 
developed a new sampling strategy that intensively selects sample units in overlapping 
areas among images acquired from adjacent rows, and automatically adds spectrally 
dissimilarity sample units from non-overlapping areas, to improve the efficiency of 
representative sample selection at the ecoregion scale. The mean OA for all scenes was 
82%. The satellite derived multi-temporal landscape quantification results revealed that 
MPB accounted for 70% of the total area of disturbance. I found that whether fire and 
MPB are linked disturbances depended on their occurring sequences. Fire severity was 
largely unrelated to pre-fire MPB outbreak severity, whereas post-fire beetle severity was 
shown to decrease with fire severity. The recovery rate varied among different 
disturbance types. Half of the clearcut and fire areas were at various stages of recovery, 
but the regeneration rate was much slower at MPB disturbed sites. Beetle outbreaks and 
fire created a positive compound effect on the seedling reestablishment, which suggests 
that beetle-killed serotinous lodgepole pines might have a new forest resilience 
mechanism to subsequent wildfire. 

Following the depiction of the disturbance pattern in landscapes, I further assessed the 
effects of a variety of biotic and abiotic factors on the outbreak dynamics in Grand 
County, Colorado. Thirty-four variables were included to develop a number of general 
linear models (GLM). Case and control samples were extracted from maps derived from 
satellite image. I first removed non-significant predictors based on the Bayesian 
Information Criterion in a multiple backward stepwise selection, and then built the model 
using the retained variables. A correction factor was added into the traditional GLM to 
account for model bias introduced by different ratios of case and control observations in 
the sample and in the population. Finally, I evaluated the model performance with an 
independent validation dataset, and generated predictive maps of MPB mortality. The 
final model had an average area under the curve value of 0.72 in predicting the annual 
area of new mortality. The results showed that neighborhood mortality, winter mean 
temperature anomaly, and residential housing density were positively associated with 
MPB mortality, whereas summer precipitation was negatively related. The extent of MPB 
mortality will expand under both RCP 4.5 and 8.5 climate-change scenarios, which 
implies that the impacts of MPB outbreaks on vegetation composition and structure, and 
ecosystem functioning are likely to increase in the future. 

Disturbance is the main driver for the heterogeneous landscape mosaic, and the 
understanding about its pattern, regime characteristics, impacts on forest resilience 
system and future trend is of great importance to many fields of research, such as carbon 
cycling, biological conservation, and environmental protection. The overall working 
approach in this dissertation provides feasible algorithms that can be applied to other 
regions, and can aid in generating consistent and high temporal frequency data on insect 
mortality and other disturbances impacting a variety of ecosystem services. 
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Chapter 1 Introduction 

1 Background 

This dissertation was motivated by the urgent need for information on current and 
projected forest disturbance dynamics to enhance and update the understanding of 
baseline carbon stocks, carbon sequestration, as well as greenhouse gas emissions, and to 
predict future potential carbon sequestration given climate, policy and land management 
scenarios in North America (Zhii et ah, 2010). Forest ecosystems in North America are 
generally considered to be a carbon sink in the past decades as a result of a faster forest 
regrowth rate than harvest rate (CCSP, 2007; Goodale et ah, 2002; Potter et ah, 2007). 
Based on estimates in forest inventories, forest ecosystems in the conterminous United 
States contain approximately 54.6 billion metric tons of organic carbon above and below 
the ground for the baseline year 1992, with an increase of 11.3 billion metric tons from 
1952 (Birdsey & Heath, 1995; Goodale et ah, 2002). Though the fluctuation of US 
carbon storage is generally considered small as reviewed from the historical records, and 
the majority of the loss is ascribed to human-induced impacts, such as forest clearing for 
urban development (Birdsey, 1992), the effects of natural disturbances should not be 
overlooked in future carbon estimation, especially in the face of climate change. 
Variation in temperature and precipitation can significantly influence the occurrence, 
timing, frequency, duration, extent and intensity of disturbances (Baker, 1995; Turner et 
ah, 1998), and ultimately alter the composition, structure and function of forest 
ecosystem (Hicke et ah, 2011; Zhou et ah, 2013). Within the United States, natural 
disturbances having the greatest effects on forests include fire, drought, windstorms, 
floodwater and ice storm, insects and disease (Barnes et ah, 1997). 

At the current stage, wildland fires receive the most research attention among all natural 
disturbance types because of their large magnitude, high severity, and their rapid 
responses to changes in climate (Dale et ah, 2001), that consequently have been 
intensively mapped and well incorporated into the disturbance modeling and emission 
estimation. While insects are another critical forest disturbance agents, their impacts on 
vegetation mortality and carbon uptake have not been well understood and considered in 
a modeling framework to date (Kurz et ah, 2008). In recent years, mountain pine beetles 
(Dendroctonus ponderosae; MPB) have caused extensive outbreaks in western North 
America, from British Columba in Canada to New Mexico in the USA, and affected 
millions of hectares of forests (Woulder et ah, 2006). It was estimated that beetles killed 
trees contributed almost equivalent carbon to those damaged by fire in the western US 
from 1984 to 2010 (Hicke et ah, 2013), which together represent 9% of the total tree 
carbon in western forests. Besides carbon, beetle infestation can affect potential fire 
hazard by altering the arrangement, composition, moisture content of forest fuels, and 
microclimate over time (Hicke et ah, 2012b; Jenkins et ah, 2008; Parker et ah, 2006; 
Schoennagel et ah, 2012). Beetle outbreaks can cause significant effects on the 
biodiversity of conifer-dominated forests by changing the nesting, roosting and foraging 
suitability of species (Martin et ah, 2006). A substantial shift in soil and water nutrient 
chemistry was observed in response to the MPB epidemics, which may degrade the 
drinking-water quality (Clow et ah, 2011; Edburg et ah, 2012; Mikkelson et ah, 2013). 
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The associated forest dieback due to beetle infestation could also impact 
evapotranspiration and albedo, alter the local surface energy balance, and finally affect 
regional temperature and climate (Boon, 2009; Maness et ah, 2012). In terms of 
economic benefits, a decrease in economic surplus is expected as a result from reduced 
timber salvage and depressed outdoor recreation industry (Abbott et ah, 2008; Prestemon 
et ah, 2013). 

Therefore, the delivery of timely and accurate information of bark beetle dynamics is 
critical in updating knowledge about ecosystem type and structure, as well as fuel 
conditions so that better model perfonnance on forest disturbance and more adequate 
assessment of forest carbon stocks and sequestration dynamics can be achieved. These in 
turn will help support climate change adaptation and mitigation activities. Such 
information is also helpful in the understanding of the disturbance regimes, their 
interactions with other disturbance types, and the effects on forest resilience. 

Besides the monitoring of contemporary beetle outbreaks, the prediction of future 
geographical distribution and range shift of this pine killer is also an important 
component in modeling carbon fluxes under various future climate change scenarios, and 
this work can hardly be accomplished without a thorough investigation of the various 
biotic and abiotic factors that govern the outbreak dynamics of bark beetles (Raffa et ah, 
2008). Although the mechanism of its population dynamics at stand or plot levels has 
been well understood (e.g., Raffa & Berryman, 1983; Logan et ah, 1998), few studies 
have focused on the investigation of landscape-level drivers that explain the spatial- 
temporal patterns of MPB mortality due to their complex interactions with fire and the 
host species as a whole in response to climate change. 

2 Review of concepts and literature 

In this section, I will first review some important concepts regarding the biology and 
population development of MPB. I will then summarize the existing knowledge and 
various techniques in detecting the beetle disturbances. Specifically, the pros and cons of 
remote sensing in forest infestation monitoring will be discussed in terms of payload 
selection, and classification methods. Finally, the interactions between insect outbreak 
and various environmental/anthropogenic factors, especially the effects of climate change 
on beetle outbreaks will be reviewed. 

2.1 Biology and epidemiology of mountain pine beetle 

As a native species to North America, MPB is a major regulator of pine forest ecosystems 
(Haack & Byler, 1993), of which the most common hosts are lodgepole pine (Pinus 
contorta), ponderosa pine (Pinus ponderosa) and western white pine (Pinus monticola), 
while lodgepole pines occupy more than 80% of the stem density of the infested stands 
(Edburg et ah, 2012). Its present geographic distribution extends from northern Mexico 
through the Pacific Northwest to northwestern Alberta in Canada, British Columbia and 
from the Pacific Coast east to the Black Hills of South Dakota, most frequently occurring 
at elevations of 1500-2600 m above mean sea level depending on the latitude (Amman, 
1973; Sambaraju et ah, 2012). 
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MPB is a univoltine species with a relatively well-defined life cycle from: egg, larva, 
pupa to adult. Adult beetles emerge and disperse between June and August, though the 
exact timing is site-specific. For each individual host tree, the infestation process starts 
from the green stage, when MPBs cut off the nutrition transport from roots by 
establishing their breeding galleries in the tree phloem, and the associated blue stain fungi 
accelerate the process by causing water stress. At this stage, the sapwood moisture will 
gradually drop but the foliage of trees remains visibly unchanged. After the death of the 
tree, the top of the crown begins to fade from green to greenish-yellow, and later to a 
uniform yellow, bright red and to brown after twelve months following mass beetle attack. 
This stage is called red stage. Normally within 3-5 years, most trees will lose all needles 
and we call that the gray stage (Amman, 1982). 

There are also evident symptoms of bark beetle attack at the stand level. In the beginning, 
the infested trees are sporadically distributed with few red or grey spots in the green stand. 
As the beetle populations grow, groups centered on the original infested trees are fonned, 
and some groups may even coalesce to large patches. At the peak of outbreaks, basically 
the whole stand appears to be effected from a landscape view. 

In the population cycle of MPB, four phases were recognized in tenns of population size 
relative to the abundance of available hosts: endemic, incipient-epidemic, epidemic and 
post-epidemic (Safranyik et ah, 2007). The beetle density in an endemic population is 
very low so that only suppressed trees with lower water and nutritional reserves can be 
attacked (Amman, 1984). Beetles at the incipient-epidemic stage are able to successfully 
initiate mass attack to a single large-diameter tree within a stand, though the clumps of 
infested trees are scattered and confined to parts of individual stands (Safranyik et ah, 
2007). Under favorable conditions, the populations will grow to epidemic level, and the 
outbreaks start taking place at the landscape scale via long-distance dispersal (Borden 
1982). In western Canada, the average duration of epidemics is around 10 years, while 
the longest recorded one continued for 18 years (Safranyik et ah, 1974). Afterwards, the 
population will decline to a post-epidemic level, as a consequence of local depletion of 
large-diameter susceptible hosts or lethal low temperatures (Safranyik et al., 1999). 

2.2 Beetle monitoring techniques 

A variety of beetle outbreak monitoring approaches exist to meet different research scales, 
extent, resolution needs, as well as management objectives. In this section, I will present 
three most common and useful streams of techniques, and discuss their benefits and 
limitations. 

2.2.1 Field survey 

Ground survey is a widely applied method to confirm disturbance agents, evaluate timber 
killed, locate trees at green stage, and collect data of the surrounding environment at the 
individual tree scale. This approach can quickly provide managers and researchers the 
current knowledge of MPB behavior at sub-outbreak levels, and has been regarded as the 
most reliable source of information about the agent responsible for forest disturbances 
(Wulder et al., 2006a). Below, I summarized three conditions when the field surveys 
behave more powerfully than other monitoring techniques: 
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1) when precise information on the physiological parameters of host trees (e.g., 
sapwood moisture, DBH), stand metrics (e.g., number of infested trees), or 
environmental variables of a particular site (e.g., soil moisture, ground fuel load) is 
required. 

2) when infestations are at the green stage, and the healthy status and beetle loads of 
individual trees should be determined. 

3) when references are needed for coarser resolution datasets, such as for aerial or 
satellite based maps (Backsen & Brian, 2013). 

However, the high economic and labor costs for field surveys restrict their application at 
a broader scale or on a more routine schedule. Because of the high per-hectare cost, field 
surveys are typically conducted on a sample basis to derive the population infonnation, 
which requires additional work on the selection of an optimal sampling strategy to avoid 
estimation bias. In addition, the long temporal gaps between plot measurements make the 
monitoring of forest disturbance status and trends challenging (Schroeder et ah, 2014). 

For instance, the revisiting cycle of a full sample of the forested landscape in the western 
US conducted by the Forest Inventory and Analysis is 10 years, which does not match the 
forest disturbance and recovery cycle in this region. 

2.2.2 Aerial survey 

Aerial survey is another method for observing forest insects and disease disturbance from 
an aircraft, delineating the affected area onto a map, and recording the associated 
attributes (such as damage type, causal agent, host, symptom). The traditional aerial 
sketchmapping is conducted using paper maps, and the advent of geospatial techniques 
enables the compatibility of aerial observations with a GIS database. The Forest Health 
Monitoring Aerial Detection Surveys (ADS) of U.S. Department of Agriculture were 
conducted as the major state and federal effort in identifying forest health status since the 
mid-twentieth century (Man, 2010), while in Canada, the Canadian Forest Service have 
conducted the Forest Insect and Disease Survey for over fifty years (McConnell et ah, 
2000 ). 

Generally, the aerial surveys have been praised for its efficiency, low cost, fine resolution, 
long temporal records and good timing in detecting insect activity. The map scale for an 
overview survey is generally 1:100,000, 1:126,720, or 1:250,000. Depending on the map 
availability, desired precision and detail, the typical coverage for an aerial survey is 
around 480 square miles per hour (McConnell et ah, 2000). The long history of archived 
sketchmapping of forest disturbance in North America promotes studies on the 
disturbance mechanism over time. Moreover, most aerial surveys are undertaken at an 
annual cycle from early-July through August, though the exact time may change 
according to different regions. This coincides with the optimum damage symptom 
expression of MPB damage as well as its univoltine feature, and thus makes it a perfect 
data source in beetle activity monitoring. 

Despite of all those merits, aerial survey data has been criticized for three main defects. 
One is the spatial imprecision at fine scales due to reasons, such as misregistration as well 
as operators’ subjectiveness (Meigs et ah, 2011), and it is most reliable when used at 
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coarse scales. A pilot study initiated in the Rocky Mountain Region demonstrated that the 
overall accuracy, producer’s and user’s accuracies as estimated from the ground truth 
points for MPB detection in lodgepole pines exceeded 80% when the spatial tolerance 
was increased to 500 m (Hohnson & Ross, 2005). The second one is the considerable 
temporal shift because of the time lag between the attack onset and the spectral signature 
changes of canopy (Meigs et al., 2011; Kautz, 2014). The third point is the perception on 
the attribute estimations. Some attributes, such as the mortality area and damage degree 
are hard to measure even during ground surveys. Thus, the sketchmapping results are not 
recommended to be extrapolated beyond reasonable expectations without more detailed 
information (British Columbia Ministry of Forests, Forest Practices Branch & Canadian 
Forest Service, Forest Health Network for the Resources Inventory Committee, 2000). 

2.2.3 Satellite-based remote sensing 

Satellite-based remote sensing platforms provide users with a wide range of choices in 
the spectral, spatial and temporal scales, and thus have been employed to address a 
variety of monitoring and decision-making issues. As a novel technique, reports on the 
mapping of bark beetle infestation using remotely sensed data were sporadic before 2000 
(e.g., Keane et al., 1994), but the growing momentum is impressive, especially in the 
most recent decade. In this part, I will focus on the comparison of beetle activity 
detection performance when different payload selection and classification methods are 
applied. 

2.2.3.1 Payload selection 

High spatial resolution imagery with its resolution less than or equal to 5 m x 5m 
(Franklin & Wulder, 2002), shows great potential for surveying in incipient population 
conditions, since this resolution of imagery captures characteristics at the stand or tree 
level (Wulder et al., 2006b). Dennison et al. (2010) used 0.5 m pan-sharpened GeoEye-1 
image to classify areas of green, red and gray canopy cover. Coops et al. (2006a) 
employed the Red-Green Index derived from 2.5 m multi-spectral QuickBird imagery to 
separate red attack crowns from the non-attack ones. An unsupervised clustering running 
on 4 m multispectral IKNONS imagery can achieve 71 %-92% accuracy in detecting low 
to medium attack (White et al., 2005). Additionally the level of detail afforded by high 
spatial resolution imagery provides critical calibration and validation data for lower 
spatial resolution remotely sensed imagery (Wulder et al., 2006b). However, because 
most high resolution payloads are operated by commercial companies, the cost per 
hectare is relatively high and the image availability is limited. Moreover, the trade-off for 
higher resolution is the smaller overall image extent, which will greatly increase the 
computational complexity and time if a broad coverage of forest needs is required. 

The spatial coverage and temporal continuity provided by coarse spatial resolution 
payloads make them suitable for monitoring forest dynamics at inter-annual to decadal 
time scales and from regional to global scales. Images with pixel size larger than 30 m 
are more frequently used in depicting general forest ecosystem change (eg. Mildrexler et 
al., 2009; Sulla-Menashe et al., 2013) than type-specific mapping, due to the less 
representation of each land cover type in a relative big image cell. But there are still some 
successful cases that take advantage of the time series trajectory. For instance, Moderate 
Resolution Imaging Spectroradiometer (MODIS) images with an 8 day revisiting cycle 
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and 500 m spatial resolution were used to map near-real-time beetle infestation in North 
America (Anees & Aryal, 2014) and gypsy moth outbreaks (Spruces et ah, 2011) through 
the use of time series. 

Because of the low availability and computational complexity of fine spatial resolution 
data, and the relative low landscape representation in coarse resolution images, moderate 
spatial resolution image is a good choice in balancing the tradeoff between accuracy and 
landscape detail. Consequently, a rising amount of studies are devoted in insect 
disturbance mapping using moderate resolution image, among which, Landsat Thematic 
Mapper (TM) data is especially popular because of its free availability, multispectral 
properties, broad spatial extent and temporal continuity. Single date Landsat image 
analysis has been successfully applied to map MPB infestations (Franklin et ah, 2003), 
but higher accuracy can be achieved when time series stacks were used (e.g. Goodwin et 
ah, 2010; Liang et al., 2014a; Meigs et al., 2011; Meddens et al., 2013; Neigh et ah, 2014; 
Skakun et al., 2003; Walter & Platt, 2013). The details about the effects of different 
classification methods will be discussed later. 

Other new generation satellite platforms, such as hyperspectral, multi-angular, light 
detection and ranging (Lidar), and radio detection and ranging (Radar) systems offer 
potential for better detection of beetle outbreaks. For example, EO-1 high resolution 
hyperspectral Hyperion moisture indices incorporating both the shortwave infrared and 
near infrared bands were significantly correlated with levels of beetle damage (White et 
al., 2007). Canopy directional reflectance obtained from the automated multi-angular 
spectroradiometer platform (AMSPEC) data was highly related with field measured 
disturbance levels (Hilker et al., 2009). 

Radar is less widely used in monitoring infestation severity or coverage, but it has shown 
promising results in studying the mass aerial migration pattern of insects through the 
boundary layer clear-air returns (Chapman et al., 2005; Nieminen et al., 2000). For 
instance, the Prince George Doppler Weather Radar was applied to estimate the aerial 
movement distance and the vertical distribution of MPB in British Columbia (Ainslie et 
al., 2011; Jackson et al., 2008). 

Lidar utilizes the round-trip times of laser pulses from signal emission to return from 
objects to measure the range between the sensor and the object surface. In forest 
ecosystem, it is mostly widely used in deriving infonnation about canopy height, density, 
and forest attributes such as biomass with higher accuracy than multispectral imagery, 
due to its greater sensitivity to forest structure variability (Hudak, 2006; Latifi et al., 

2010). The typical function of lidar in beetle-related research is to relate lidar metrics 
with field observations (Bright et al., 2013), or with disturbance history metrics derived 
from multi-spectral images (Bright et al., 2012; Pflugmacher et al., 2014), or with high 
resolution digital camera imagery (Bater et al., 2010; Coops et al., 2009) to map percent 
dead basal area, aboveground biomass or carbon stocks. Moreover, fine-scale lidar data 
can be used as a reference to validate the results of coarser resolution image 
classifications (Bright et al., 2014). 


6 



In summary, factors to be considered in the selection of payloads when planning to map 
MPB infestations include the spatial resolution, timing, frequency, extent, cost, 
computational complexity, thematic detail and accuracy. It is virtually impossible for any 
single remote sensing platform to fulfill all these requirements. Users need to rank the 
factor importance and well understand the payload characteristics in order to match the 
information needs of forest managers, decision makers or research in other scientific 
domains. 

2.2.3.2 Classification methods 

A variety of classification methods have been proposed to map tree mortality. According 
to the classification units, the methods can be grouped into pixel-based and object- 
oriented approach (eg. Bunting & Lucas, 2006; Coggins et ah, 2008). Depending on the 
density of available images, they can be categorized into single-date and multi-temporal 
classification. 

Single-date classification is the most applied technique, because of its easy operation and 
shorter preparation time needed to conduct nonnalization and image registration. This 
approach is suitable for urgent updates on the current infestation condition. Nonetheless, 
two major gaps exist in terms of long-term tree mortality detection. First, most current 
single-date classifications are developed using high resolution images, such as 30 cm 
resolution aerial imagery (Meddens et ah, 2011), 0.5 m resolution GeoEye-1 (Dennison et 
al., 2010), 2.5 m resolution Quickbird (Coops et al., 2006a), and 4 m IKONOS (White et 
ah, 2004). Although technically feasible, it is not feasible to obtain fine resolution images 
consistently over time. Second, traditional single-date machine learning algorithms 
require large training sample size to represent all the variability present in a category. 
However, the acquisition of adequate time-series ground truth survey data is usually 
impossible considering the high economic cost and time efficiency. Some studies would 
gather training and test samples based on the images to be classified, but this way has 
been criticized for introducing substantial uncertainties (Powell et al., 2004). 

Multi-temporal classification has been demonstrated to be more effective in separating 
multiple post-outbreak attack stages from other classes than single-date classification 
(Byme et al., 1980; Meddens et al., 2011). An integrated approach using multi-temporal 
remotely sensed data has begun to emerge (Woulder et al., 2006; Hilker et al., 2009; 
Coops et al., 2010; Meigs et al., 2011). In general, the multi-temporal classification 
methods can be categorized into three main types: threshold-based, trajectory-based, and 
post-classification change detection. 

1) Threshold-based change detection employs certain thresholds that represent 
different degrees of crown foliage changes to separate the mortality from no change 
in either manual or machine learning algorithms (eg. Wulder et al., 2006b). Under 
most conditions, indices transformed from the image digital number or spectral 
reflectance were used for general surveys of beetle damage. Some typical indices 
include: Red-green index - the ratio of red to green reflectance (Coops et al., 2006a; 
Hicke & Logan, 2009); Enhanced wetness difference index - subtraction of the 
wetness component in the Tasseled Cap Transformation (Cohen et al., 2003; 

Franklin et al., 2003; Skakun et al., 2003). Threshold-based change detection can be 
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easily operated and interpreted when displaying the difference in wetness in colors. 
However, the cut-off values are site-specific since they highly depend on the 
landscape characteristics, such as the infestation status, heterogeneity, tree cover and 
topography. A uniform value could lead to potential bias if it has been applied across 
a relative large extent. In addition, while stand-replacing events are featured by 
prominent changes in vegetation structure and function that can be easily detected 
with a large cut-off value (Skole & Tucker, 1993), relatively subtle changes arising 
from insect attacks might be harder to separate from background changes using a 
small threshold. 

2) Trajectory-based mapping detects changes using the temporal trajectory of 
spectral values through simultaneous analysis of all available images in the dense 
image stack. Before the public open access to the Landsat archives, most studies 
employed the high temporal frequency but coarse resolution satellite data such as 
AVHRR and MODIS (Myneni et ah, 1997; Coppin et ah, 2004; Lunetta et ah, 2006). 
Now, a rising number of studies are working on the algorithm development of 
Landsat time-series stack in monitoring forest disturbance history. Two successful 
examples include the Landsat-based detection of trends in disturbance and recovery 
(LandTrendr; Kennedy et ah, 2010) and the Vegetation Change Tracker (Huang et 
ah, 2010). The core idea of LandTrendr is temporal segmentation, that decomposes 
often noisy time series into a sequence of straight-line segments with distinct 
features designed to capture the event that happened for the duration of each segment. 
VCT works similar to LandTrendr with differences in that each individual image is 
first analyzed to create a mask for water, cloud and cloud shadows, as well as forest 
samples, and the time-series of derived indices and masks is then used to 
characterize changes (Huang et ah, 2010). Both methods recognize disturbance 
events by the rate of changes in time series of spectral properties. This makes the 
head and tail years of the trajectory more difficult to judge, because the attribution 
for a particular year is entirely detennined with regard to its spectral deviation from 
the two adjacent years (Kennedy et ah, 2010). In addition, the event labeling 
depends on pre-implemented modules in the package, and thus is not flexible in 
adjusting to user-defined functions. Furthermore, both algorithms are able to find 
most of the abrupt disturbances, but neither was fully optimized to capture the 
chronic and low severity insect damage (Schroeder et ah, 2014). Nonetheless, this 
type of method is promising because of the free availability of more satellite images 
at fine to medium spatial resolution in the future, and temporally dense observations 
are necessary to distinguish between type and frequency in heterogeneous 
landscapes (Neigh et ah, 2014). 

3) Post-classification change detection: Unlike the fonner two types of methods that 
focus on the relative differences in spectral response between images covering the 
same extent or the changing trend of the spectral curve over time, post-classification 
mapping runs classifiers on each individual image first, and then compares the 
results to detect changes in land cover types. This type of method has been widely 
applied in various land cover change or forest disturbance studies (eg. Kuemerle et 
ah, 2009; Griffiths et ah, 2013). Besides the final change products, it also generates 
the land cover map for each image, which could be attractive to managers or 
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decision makers if they need consistently updated land cover synthesis. However, 
because the accuracy of change detection highly depends on the classification results 
at each date (Coppin et ah, 2004), the errors could propagate and lead to poor 
accuracy, especially under the circumstance when longer time series is involved. 
Some methods were proposed to alleviate this issue and the common principle lies in 
the development of a land cover transition matrix, since the consequence of error 
propagation is the occurrence of illogical land cover transitions. Liang et al. (2010) 
identified the uncertainties in the MODIS Collection 5 global land cover products 
using a land cover type conversion lookup table. Liu and Cai (2012) employed the 
spatial-temporal contextual information modeled by Markov random field theory to 
reconstruct the change trajectories. Wang et al. (2014) applied the similar approach 
to carry out the global land cover mapping at 250 m resolution with multiple year 
time series MODIS data. Those efforts greatly improved the accuracy of post¬ 
classification accuracy, but the design of the land cover transition matrix is mostly 
arbitrary and the intricacies of the algorithm as well as the computational complexity 
might limit its application outside of the remote sensing community. 

2.3 Climate change and mountain pine beetle 

Bark beetles are cold blooded. They can respond quickly to changing climate by shifting 
their geographical distribution to maintain their climatic niche (Carroll et al., 2003), 
which makes climate warming the biggest concern in affecting the extent and severity of 
beetle disturbance (Bale et al., 2002). Some direct effects of climate change on bark 
beetle include developmental timing and cold tolerance, and indirect effects include 
community associates, host-tree physiology, and host-tree distribution (Bentz 2010). 
Evidences of range expansion into new habitats as a result of increasing temperature have 
been documented in many cases (eg. Battistia et al., 2006; Nealis & Peter, 2009), and 
they are attacking hosts farther north and at higher elevations than historic nonns (Tran et 
al., 2007; Safranyik & Carroll, 2006). 

The bark beetle’s life cycle is temperature dependent (Figure 1). Their daily emergence 
rates are proportional to cumulative degree-days above 14.4 °C, and captures of released 
MPB were directly related to heat accumulation above 16 °C (Safranyik et al., 1989). 

Most beetles fly between the temperature from 22°C to 32 °C (Safranyik, 1978), and the 
flight propensity increases with ascending light intensity and relative humidity if 
temperature is within the optimal range. Once the temperature exceeds 35 °C, the flight 
propensity begins to decline and will disappear above 38 °C (McCambridge, 1971). Cold 
temperature is the most threating factor of beetle mortality (Safranyik, 1978; Cole, 1981): 
eggs can not survive below -18 °C (Reid & Gates, 1970), and the lethal temperature for 
pupae is between -18 °C and -34 °C, and for adults is between -23 °C and -34 °C (Logan et 
al., 1995). The degree of cold tolerance is determined by the amount of glycerol (Bentz & 
Mullins, 1999), which gets accumulated with the mature of larvae, and usually reaches its 
peak during the winter period from December to February. Thus, if an extreme cold 
weather that lasts for several days strikes before the sufficient accumulation of glycerol 
or after the beetles have already metabolized it, massive mortality can occur (Safranyik et 
al., 1974). 
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Although some studies suggest that there may be sufficient genetic variability in MPB 
populations to match changes in the climatic environment (Bentz et al., 2001), the rate of 
climatic niche evolution is not predicted to be faster than that of climate change (Quintero 
& Wiens, 2013). This is the basis for many studies in estimating the potential distribution 
of bark beetle under current and projected future climate conditions. Both process-based 
and statistical models have been applied in exploring the responses of beetles to changing 
climate and estimating the future potential habitats. Process-based models are focused on 
linking population dynamics with climate in a physiologically explicit way, whereas 
statistical models build the association of outbreaks with climate events from the 
perspective of data correlation. 

Bentz et al. (2010) used BioSIM that was driven by hourly temperature to identify 
regions with a high potential for MPB outbreaks across North America. Using a process- 
based model, Hicke et al. (2006) found that projected wanning in the western United 
States will result in substantial reductions in the overall area of adaptive seasonality. 
Aided by an ecological niche model, Evangelista et al. (2011) predicted that new areas of 
forest susceptible to MPB mortality would emerge over time but the existing distribution 
would also shrink under future climate condition, leading to an overall decrease in 
suitable habitats. Based on the results of logistic regression models, Sambaraju et al. 
(2012) indicated that increases in mean temperature by ID to 40 would profoundly 
increase the outbreak risk at higher elevations first then at higher latitudes, whereas 
variance in mean temperature would not change the overall trend in outbreak potential. 
Carroll et al. (2004) modeled the climatic suitability using historical weather in British 
Columbia, Canada, and a substantial shift in climatically suitable habitats was found 
towards northern and higher elevations. Noticeably, the MPB outbreaks from historical 
survey data followed this shifting pattern. In Oregon and Washington, USA, the odds of 
outbreak at locations with highly suitable weather was 2.5 times higher than those with 
low suitability, as calculated from a weather suitability index (Preisler et al., 2012). 
Although a large body of studies have confirmed the importance of climatic variables in 
shaping the outbreak dynamic patterns at the macro-scale (eg. Chapman et al., 2012; 
Jewett et al., 2011; Lester & Irwin, 2012; Mitton & Ferrenberg, 2012), many of them 
only consider a limited number of explanatory variables because of model complexity, 
especially for the process-based approaches, which might eventually introduce errors into 
the simulated model results. 
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Temperature (°C) 

Figure 1. The relationship between emergence frequency of mature MPB and 
temperature (adapted from McCambridge 1971). 

3 Objectives and thesis structure 

As revealed from the literature review in Section 2, although great improvements have 
been achieved in terms of MPB monitoring, detection and modeling, caveats and 
knowledge gaps still exist. The overarching goal of my dissertation is to enhance the 
time-series mapping ability of medium resolution remotely sensed data in areas 
characterized by heterogeneous disturbance-derived mosaic and rapid landscape pattern 
changes. This will allow for a better and more seamless integration with ecological 
modeling as well as landscape analysis, to improve our understanding about the resilience 
mechanisms of forest ecosystem and facilitate the forest management and decision 
making. Each chapter addresses one specific objective, and the structure is organized as 
follows: 

Chapter 2 developed a practical workflow to derive the progressive MPB outbreaks in a 
time sequence using medium resolution remotely sensed data over a decadal long period. 
I proposed a forest growth trend analysis to conduct the mapping by eliminating the 
observational noise in the time-series, capturing the disturbance events as represented by 
spectral segments, and quantitatively attributing events based on their distinct features, 
with an integrated Landsat temporal trajectories and decision tree techniques. The 
accuracy was compared with two single-date classifiers, i.e., the Maximum Likelihood 
classifier and Random Forest with varying number of training samples, to assess the 
sensitivity of classification outcomes to the sample size. 

Chapter 3 advances the procedure in Chapter 2 and conducts the long-term mapping in 
the Southern Rocky Mountains ecoregion, with the aim of depicting the disturbance and 
recovery patterns, characterizing their regimes, and describing the disturbance- 
succession pathways. I first developed a sampling strategy that took advantage of the 
Landsat scene-overlapping area between images to promote sampling efficiency at a 
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large scale, and then automatically added representative samples based on the spectral 
similarity rule to improve the mapping accuracy. The analysis method proposed in 
Chapter 2 was improved via an automatic event attribution step using random forest and a 
set of ecological meaningful features that contain the spatio-temporal information of the 
focal segment. Based on the generated annual time-series of disturbance/recovery, I 
derived their key regime parameters, such as onset year, magnitude, duration, and 
disturbance-recovery pathways, to investigate the complex interactions among various 
disturbance types in this region, and the effects on the forest resilience mechanism. 

In chapter 4,1 quantify the anthropogenic and environmental drivers that explain the 
dynamic patterns of MPB mortality at the landscape level, and simulate areas with future 
potential MPB mortality under projected climate-change scenarios. I applied the general 
linear models (GLM) to rank the importance and significance of various biotic and 
abiotic factors. The response variables were defined as the new MPB mortality, and the 
case and control samples were extracted from the remote sensing derived current 
distribution of pine beetle outbreaks in a stratified random sampling manner, of which the 
sample sizes were determined via Moran’s I to avoid spatial autocorrelation effects. I 
then complied a set of comprehensive landscape-level explanatory variables, including 
anthropogenic, topographic, vegetation, climate and dispersal-related parameters from 
various geospatial sources. The GLM with the best performance, which was evaluated in 
both qualitative and quantitative ways, was applied to predict potential new areas of MPB 
mortality using the outputs of downscaled global climate models from the Coupled 
Model Intercomparison Project 5. 

Finally, in chapter 5,1 provide conclusions from this research and point out some future 
research directions that can be carried out based upon the outcomes of this dissertation. 
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Chapter 2 Trajectory-based classification for beetle outbreak detection 
with time-series Landsat stack 

This article has been published previously and is reproduced here with permission from 
the publisher, Elsevier 

Liang L, Chen YL, Hawbaker T, Zhu ZL, Gong P. 2014. Mapping mountain pine beetle 
infestation through growth trend analysis of time-series Landsat data. Remote Sensing, 
6(6): 5696-5716. 


Abstract 

Disturbances are key processes in the carbon cycle of forests and other ecosystems. In 
recent decades, mountain pine beetle (MPB; Dendroctonus ponderosae) outbreaks have 
become more frequent and extensive in western North America. Remote sensing has the 
ability to fill the data gaps of long-tenn infestation monitoring, but the elimination of 
observational noise and attributing changes quantitatively are two main challenges in its 
effective application. Here, we present a forest growth trend analysis method that 
integrates Landsat temporal trajectories and decision tree techniques to derive annual 
forest disturbance maps over an 11-year period. The temporal trajectory component 
successfully captures the disturbance events as represented by spectral segments, whereas 
decision tree modeling efficiently recognizes and attributes events based upon the 
characteristics of the segments. Validated against a point set sampled across a gradient of 
MPB mortality, 86.74% to 94.00% overall accuracy was achieved with small variability 
in accuracy among years. In contrast, the overall accuracies of single-date classifications 
ranged from 37.20% to 75.20% and only become comparable with our approach when 
the training sample size was increased at least four-fold. This demonstrates that the 
advantages of this time series work flow exist in its small training sample size 
requirement. The easily understandable, interpretable and modifiable characteristics of 
our approach suggest that it could be applicable to other ecoregions. 

1 Introduction 

Understanding patterns of forest disturbances is important for assessing past and future 
ecosystem structure, function, productivity and diversity (Turner, 2010; Wang et al., 
2012). In recent years, mountain pine beetle (MPB) outbreaks have occurred over 
extensive areas, from British Columba in Canada to New Mexico in the United States. 
MPB is endemic to North America; however, the recent outbreaks have reached epidemic 
levels that have affected millions of hectares of forests (Hieke et al., 2013). The effects of 
MPB infestation range from altered surface fuel and wildfire hazards (Hicke et al., 
2012b; Jenkins et al., 2008; Parker et al., 2006; Schoennagel et al., 2012), changed 
vegetative composition (Collins et al., 2011), converted live carbon sinks to dead and 
slowly decaying carbon sources (Caldwell et al., 2013; Kurz et al., 2008; Running, 2008), 
impacted nutrient cycling and water quality (Edburg et al., 2012; Mikkelson et al., 2013), 
shifted evapotranspiration and albedo, modified local surface energy balance (Boon, 
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2009) and changed regional climate (Maness et al., 2013). Those effects are predicted to 
increase as a consequence of the direct and indirect effects of climatic changes (Bentz et 
ah, 2010; Raffa et ah, 2008; Williams et ah, 2010). However, unlike wildfires that have 
received intensive mapping efforts at different scales, characterizations of the extent and 
severity of MPB using multi-temporal data have only begun to emerge (Coops et al., 
2006b; Goodwin et al., 2010; Meddens et al., 2011; Meigs et al., 2011; Wulder et al., 
2006a). 

Remote sensing, both with aerial photography and satellite imagery, has been widely 
acknowledged as a useful way to map tree mortality over large areas (Hansen & 
Loveland, 2012). Landsat time series stacks (LTSS) are well suited to this type of 
analysis, because they are freely available, cover a long time period, have a broad spatial 
extent, are multispectral and have temporal continuity (Franklin et al., 2003; Wulder et 
al., 2006a). Two frequently used disturbance mapping approaches are (1) trajectory- 
reconstruction and (2) decision trees (DTs). Trajectory-based change detection examples 
include the Landsat-based detection of trends in disturbance and recovery algorithm 
(LandTrendr) (Kennedy et al., 2010) and the Vegetation Change Tracker (Huang et al., 
2010; Masek et al., 2013; Vogehnann et al., 2011). Both methods recognize disturbance 
events by the rate of changes in time series of spectral properties. This makes the first and 
last years of the trajectory more difficult to judge, because the attribution for a particular 
year is entirely detennined with regard to its spectral deviation from the two adjacent 
years (Kennedy et al., 2010). In addition, the event labeling procedure depends on pre¬ 
implemented modules in the software and, thus, is not flexible in adjusting to user- 
defined functions. DT approaches include work by Goodwin et al. (Goodwin et al., 2008) 
who utilized the annual changes in normalized difference moisture index and a set of 
thresholds to classify the presence of MPB infestation over 14 years, resulting in overall 
accuracies ranging from 71% to 86%. DT-based approaches have the capacity to 
integrate expert knowledge and flexible classification schemes. However, poor image 
quality from sensor drift, geometric misregistration and topographic and cloud shadows 
will introduce distortions to images and affect the accuracy of change detection. 

In this study, we describe a forest growth trend analysis method that integrates the robust 
change detection strengths of the trajectory reconstruction approach and the flexible user- 
defined advantages of DT approaches to characterize MPB mortality with an 11-year 
LTSS. We were interested in determining whether or not this proposed method could 
reach an accuracy target of 85% or greater. Additionally, we compared its accuracy with 
single-date image classifications, which are still important and commonly used in 
disturbance mapping, especially in two-date post-classification change detection. We are 
aware that the accuracy of single-date image classification is not only determined by the 
input variables and classifiers, but is also heavily dependent on the quality and quantity of 
training samples (Gong et al., 2013), which may have larger impacts than the implemented 
techniques (Campbell, 2002). Thus, we also assessed the effect of different sample sizes 
on the accuracy of single-date image classifications, and we tested how much our 
proposed method can improve sampling efficiency as compared with the single-date 
classifiers. 

2 Study area and data 
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2.1 Study area 

Our study area was located in Grand County, in north-central Colorado (Figure 1), with a 
total land area of 4783 km 2 . This landscape is dominated by evergreen forests mainly 
composed of lodgepole pine (Pinus contorta), along with other tree species, as well as 
shrub and grass cover types. Summer in Grand County is typically hot and dry, which is 
suitable for MPB development. The average precipitation per year is approximately 400 
mm, whereas the U.S. average is 930 mm. This county has one of the highest 
concentrations of lodgepole pines and has been heavily affected by the MPB outbreak 
that started around 2002 and reached the peak of tree mortality between 2005 and 2008 
(Klutsch et ah, 2009). 
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Figure 1 . The study area used for our analyses. The Southern Rocky Mountain ecoregion 
is shown in blue, and Grand County, Colorado, USA, is shown in red. The imagery is a 
false-color composite of one Landsat image acquired on 21 August 2009. 
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2.2 Landsat time series stack 


Grand County lies entirely within Landsat path 34 row 32. We restricted the image 
acquisition season from June to September, because this snow-free period is the 
recommended time window to detect MPB mortality (Wulder et ah, 2006a). For each 
year, one cloud-free or mildly contaminated image was selected to build the LTSS, and if 
no such image existed, several partly cloudy images with a similar phenology were used 
as replacements (Table 1). All data were processed to Level IT (terrain corrected). 
Radiometric calibration and atmospheric correction for cloud and aerosol effects were 
perfonned using the Landsat Ecosystem Disturbance Adaptive Processing System [34], 
Clouds and their shadows were automatically screened using the Function of mask 
(FMASK) algorithm (Zhu & Woodcock, 2012). Cloud removal effects were judged 
satisfactory via visual checking, though the extent of cloud cover was slightly 
overestimated. We accepted this false positive error, because commission errors were 
more desirable than omission errors in cloud filling. The number of clear observations in 
the LTSS over the period of 2000 to 2011 is shown in Figure 2, and 98% of the study 
region contained more than 10 cloud-free observations. 


Table 1. Characteristics of the Landsat images used in this analysis. 



Year 

Month and 

Julian 

Cloud Cover 

Cloud Cover over MPB Host 


Day 

Day 

(%) 

Area (%) * 

1 

2000 

Jul 13 

193 

42.35 

1.09 

2 

2000 

Jul 21 

201 

4.74 

3 

2001 

Jun 28 

179 

18.87 

4.24 

4 

2001 

Jul 22 

203 

17.05 

5 

2002 

Jul 1 

182 

1.92 

0.15 

6 

2003 

Jul 4 

185 

5.48 

0.49 

7 

2004 

Jul 8 

188 

10.49 

7.94 

8 

2005 

Sep 11 

254 

2.08 

0.25 

9 

2006 

Jul 28 

209 

2.30 

1.26 

10 

2007 

Jul 31 

212 

20.39 

4.22 

11 

2007 

Aug 16 

228 

13.19 

12 

2008 

Aug 20 

231 

28.91 

9.79 

13 

2008 

Sep 5 

247 

36.83 

14 

2009 

Aug 21 

233 

2.04 

0.00 

15 

2010 

Sep 25 

268 

3.37 

0.89 

16 

2011 

Aug 11 

223 

22.13 

2.67 

17 

2011 

Aug 27 

239 

16.45 


Notes: MPB, mountain pine beetle; * the final cloud cover in a particular year in 
MPB host forest areas is after using clean pixels in one image to fill cloudy areas in 
another image that is in the same phenological state. 


Since MPBs attack a limited number of tree species, we confined our analysis to the 
extent of their host species in the LANDFIRE Refresh 2001 Existing Vegetation Type 
data layer (EVT), which was developed using circa 2001 (1999-2003) Landsat imagery 
(Nelson et ah, 2013). The spatial resolution of the LANDFIRE data matches the LTSS, 
and the time period of the imagery from which the EVT data were derived depicts the 


16 





v.oZoor s.oioor \.o«ot \.o$«6C \.or«6( 


pre-disturbance distribution of vegetation in Grand County. Currently, MPB mortality is 
concentrated in lodgepole and ponderosa pine (Pinus ponderosa ) forests, but other pine 
species, like limber pine (Pinus flexilis) and bristlecone pine (Pinus aristata ), have also 
been found to be suitable hosts (Gibson et ah, 2008). We included all EVT pixels that 
contain any of those species (Table SI), and smoothed the image with a 3 x 3 majority 
fdter, which was designed to reduce observed salt-and-pepper effects and to minimize 
omission errors that are much higher than commission errors for the selected classes 
(National c2001 Assessment, 2001). 
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Figure 2. Forest type map for Grand County, from (a) the LANDFIRE existing 
vegetation type data and (b) the number of clear Fandsat observations over the period of 


2000 to 2011. 


2.3 Reference sample selection 

We randomly placed training samples over forested areas and reconstructed their 
disturbance history with the aid of normalized burn ratio (NBR) time series from the 
LTSS and the 1-m resolution U.S. Department of Agriculture (USDA) National 
Agricultural Imagery Program (NAIP) imagery (obtained from the USDA Geospatial 
Data Gateway), which were available for years 2005, 2009 and 2011 in Grand County. 
The NBR time series were calculated using bands 4 and 7 from the Fandsat imagery (Key 
& Benson, 2006) and then multiplied by 100, so they could be represented as integer 
values instead of floating point values. For each training location, we identified the 
disturbance event type, duration and onset time (Table 2). Disturbance event types 
included: healthy, MPB mortality or clearcuts. Fire was not included, because large fires 
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did not occur in our study area during the time period of this study (2000 to 2011), as 
indicated by the U.S. Geological Survey Monitoring Trends in Bum Severity database 
(Eidenshink et ah, 2007). In total, we collected 106 training samples, among which, 53 
were sites with MPB mortality, 37 were clearcuts and 16 were healthy. We selected the 
relatively small number of training samples, because many studies suggested that 
trajectory-based methods perform well when trained with a relatively small sample size 
(Kennedy et ah, 2010; Meigs et ah, 2011; Sulla-Menashe et ah, 2013). Statistical 
classifiers require more training data to fully represent each class in the feature space 
(Mathur & Foody, 2008). Thus, we expanded the original training sample size by 10 
times to a total of 1,060 training samples. For those newly added samples, we only 
identified their disturbance types in the three years for which NAIP imagery were 
available instead of the whole NBR trajectory from the FTSS. 

For validation purposes, we randomly sampled points in strata that were formed in the 
complete temporal-spectral space as a representation of the gradient of MPB mortality. In 
this manner, our test samples were able to cover a complete spectrum of disturbance 
events and a range of mortality. Considering the large data volume of the 11-year FTSS, 
for which the dimension was 2964 (rows) x 2901 (columns) x 102 (bands), we employed 
an efficient unsupervised classification tool, Clustering Based on Eigenspace 
Transformation, to improve the speed of clustering over conventional K-means without 
losing accuracy (Chen & Gong, 2013). Fifty strata were developed, and we 
proportionately allocated a total of 100 points to each stratum. In the absence of sufficient 
ground and ancillary information, we visually interpreted the status for each point at each 
individual year. By checking the proportion of points in the three classes, we found that 
they are almost equally distributed without severe over- or under-representation; the mean 
number of validation samples in each disturbance type each year was 39 (healthy), 38 (MPB 
mortality) and 23 (clearcuts; Table S2). 


Table 2. Examples of normalized burn ratio (NBR) trajectories, National 
Agricultural Imagery Program (NAIP) data and in situ photos for healthy forest, 
MPB mortality and clearcut classes. 

Class Landsat NAIP Photo 
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3 Method 


3.1 Temporal segmentation 

We developed time series of disturbance maps by combining the temporal trajectory 
analysis approach of LandTrendr (Kennedy et ah, 2010) with a decision tree labeling 
procedure. Temporal trajectory analysis was used to decompose often noisy time series 
into a sequence of straight-line segments with distinct features designed to capture the 
event that happened for the duration of each segment. NBR (Key & Benson, 2006) was 
chosen to represent the vegetation condition, because of its sensitivity to disturbance 
(Kennedy et ah, 2010) and its capacity to minimize radiometric differences among 
images. The general procedure was to use a set of parameters to enhance the signal-to- 
noise ratio of the NBR time series in order to capture the salient events happening in each 
period or segment. Details can be found in Kennedy et al. (2010), and the final effect of 
temporal segmentation, whether it was under-segmentation (split up into too few parts) or 
over-segmentation (subdivided into too many parts), was detennined by the parameters 
listed in Table 3. Optimal LandTrendr parameter settings have been tested in the forests 
of the Northwest Pacific Region at both Landsat and MODIS scales (Kennedy et al., 
2010; Sulla-Menashe et al., 2013). However, because we are not aware of reported 
parameter results in the literature specific to the Southern Rocky Mountain ecoregion, we 
tested a range of candidate parameter values (Table 3). 


Table 3. LandTrendr parameters, descriptions, suggested values from Kennedy et al. 
(2010) and the values tested for this study and ultimately selected for this study to reduce 
the raw normalized burn ratio (NBR) time series values to a number of linear segments. 


Parameter 


Despike 


Max_seg 


V ertexOvershoot 


pval 


Kennedy 

Description and Units et al. 

_ ( 2010 ) 

An outlier will be removed if the 
proportional difference in NBR ^ 

values between two adjacent points 
is less than this parameter. 

The maximum number of segments ^ 

allowed in fitting. 

The first round regression detected 
vertices can overshoot (max_seg + 3 

1) vertices by this value. 

A pixel will be considered as no 
change if its /.i-valuc of the F-statistic 0.05 

for the entire trajectory is above this 


Values Value 
Tested Selected 


0.75; 0.9; 

1 


0.90 


3; 4; 5 4 


0; 3 0 


0.05; 0.1 0.1 
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threshold. 

Recovery threshold 

If the slope of a segment positively 
spans the whole spectral range 
within 1/recovery threshold years, it 
will be rejected. 

0.25 

0.25; 1 

0.25 


A simpler model will be chosen if its 




BestModelProportion 

F-statistic exceeds this proportion of 
the best fitting model. 

0.75 

0.75; 1 

0.75 


We ran LandTrendr on all combinations of the candidate parameter values, which 
resulted in 144 unique combinations in total, and then visually compared their temporal 
segmentation results with the original NBR time series for the 106 training samples. Thus, 
the calibration involved 15,264 comparisons. Each individual comparison for a sample was 
coded in a binary state using two criteria: the overall match in the shape and the timing of 
vertices. A good match occurred if the disturbance event that happened at that pixel was 
successfully captured, and the slope and magnitude for each segment deviated little from 
the central trend. The vertex timing was judged by how well the algorithm captured the 
onset and ending time of each disturbance event. If both criteria were satisfied, that 
comparison would be scored as 1, and otherwise, it would be 0. The optimal parameter 
set was decided based upon the pooled rankings across all training samples. The 
processing time took about 200 minutes (on 5,376,538 unmasked pixels) for this step, 
which was run on a workstation with an Intel i7-2600 3.4 GHz quad core CPU and 8 
gigabytes of physical memory. 

3.2 Decision tree-based spectral segment labeling 

The outputs of the temporal segmentation included the fitted segments, and a set of 
vertices connecting individual segments. We employed a DT to attribute disturbance 
events associated with those vertices and segments, with leaf nodes representing class 
labels and branches representing decision scenarios. The DT was built with predictors 
characterizing each segment: the disturbance occurrence onset time, duration and the 
magnitude of NBR change. We defined the onset time as the starting year of a segment. 
Magnitude was the absolute change in NBR between the ending year and the onset year, 
whereas duration was the time elapsed over the length of one segment. Different types of 
disturbance events were characterized by distinct features. For instance, clearcuts could 
be recognized by segments with a short duration and large, negative change in NBR, 
whereas MPB mortality segments tended to have a gradual, but long-duration decline in 
NBR values. 

Decision rules were used to identify forest conditions in a top-down sequence (Figure 3), 
and the definition of key parameters are explained in Table 4. These parameters were 
manually calibrated with the same training samples described in Section 3.1. Specifically: 

(1) Segments were separated into regrowth, stable or disturbed status according to 
their magnitude (mag) at the first level of the DT. If the absolute magnitude of one 
segment was less than the pre-set threshold (<thre_magl), it was treated as 
disturbed. If the magnitude was within a threshold range (±thre_magl), then it was 
considered to be stable. Otherwise, it was considered to be in a regrowth stage 
(>thre_magl). 
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(2) To label a segment as regrowth or stable, healthy vertices were identified based on 
an NBR threshold (thre_vertexl). Thus, a vertex in either a regrowth or stable 
period with its NBR value greater than the thre_vertex 1 parameter was labeled as 
healthy. 

(3) Segments classified as disturbed in Step 2, were labeled as clearcuts if their rate of 
change in NBR was greater than a second threshold (thre_mag2); otherwise, 
segments were labeled as MPB mortality. The rate of change was defined as the 
average NBR change per year, which is a reflection of both magnitude and 
duration. The thre_mag2 parameter was designed based on the assumption that 
clearcuts always have more abrupt and rapid decreases in NBR than MPB 
mortality. 

(4) In the third level of the DT hierarchy, the label from the previous year (pre label) 
was critical in determining the label for the current vertex. This was based on the 
assumption that events are temporally dependent and forests can only logically 
transition from certain states to another. For instance, clearcuts often result in 
abrupt declines in NBR followed by persistent, but slow increases in NBR. 
Although the magnitude of change in the post-clearcut period is not as sharp as 
that of the pre-clearcut period, the vertices in the subsequent years will still be 
assigned as clearcut to ensure that the disturbance class labels have temporal 
consistency. 

(5) The final step involved attributing the vertices for the first year of the time series. 
We separated this step, because there was no vertex infonnation prior to the first 
year. After MPB mortality, standing tree tru nk s and branch residuals remain on 
site, whereas clearcuts usually have a significant amount of bare ground. Thus, 
NBR values in areas with MPB mortality should generally be higher than they are 
in clearcut areas. The thre_vertex2 parameter defined the cutoff value for 
separating clearcuts and MPB mortality in these cases. 
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Vertex 



Figure 3. The decision tree used to attribute temporal segments into the different 
disturbance classes. 


Table 4. Abbreviations and control parameters used in the decision tree. 


Parameter 

Description 

mag 

Magnitude of change in NBR between the current vertex and the following 

V ertex : NBRnext year NBRcurrent year. 

cur 

The NBR value for current vertex. 

pre label 

The vertex label in the previous year. 

thre magi 

The magnitude threshold that distinguished stable from either disturbance 

or recovery. 

thre mag2 

The magnitude threshold that separated MPB mortality from clearcut. 

thre vertex 1 

If the current vertex NBR value was above this threshold, it was treated as 
healthy, otherwise, it was further analyzed into disturbance types. 

thre vertex2 

If first year’s vertex value was below this threshold, it was treated as a 
clearcut pixel and MPB mortality otherwise. 

Notes: MPB - 

mountain pine beetle; NBR - nonnalized bum ratio. 


3.3 Post-labeling process 

Both the temporal segmentation and DT were run at the pixel level. Thus, it was possible 
that two adjacent pixels could simultaneously experience the same event, but were 
labeled differently due to errors. Since most disturbance events were patchy, the behavior 
of spatially adjacent pixels could be used to enhance the robustness of the classification. 
Based on this logic, we filled small holes in a large patch that were undetected or 
mislabeled. To do this, we first applied a 3 x 3 majority filter to reduce the salt-and- 
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pepper effect and merged pixels into more complete patches, and then, we used a 3-year 
temporal filter to remove illogical transitions over time. For example, this step would 
change a pixel classified as MPB mortality in Year 1, healthy in Year 2 and MPB 
mortality in Year 3, to MPB mortality for all 3 years. 

3.4 Single-date classification 

Single-date mapping strategies were assessed using a traditional maximum likelihood 
classifier (MLC) and a random forest (RF) algorithm. MLC is a parametric method 
widely used in forest remote sensing (Breiman, 2001). In contrast, RF is a machine- 
learning technique that combines the results of many (thousands) of weak classifiers, 
such as classification and regression trees (Dietterich, 2000). RF is often praised for its 
good overall perfonnance, as well as robustness to noise (Tucker, 1979). We 
implemented the classifiers on both the raw spectral bands and a combination of raw 
spectral bands and spectral indices that were sensitive to forest defoliation and canopy 
water content (Table 5, Crist & Cicone, 1984; Gao, 1996; Key & Benson, 2006; Johnson 
& Ross, 2005; Rock et ah, 1986; Wilson & Sader, 2002). 


Table 5. Spectral indices used in the single-date classifications, including fonnulas 
and references. 


Spectral Indices 

Formula 

Reference 

Normalized Difference 
Vegetation Index 

(b4 - b3)/(b4 + b3) 

Wilson & Sader, 2002 

Nonnalized Difference 
Moisture Index 

(b4 - b5)/(b4 + b5) 

Rock et al., 1986 

Moisture Stress Index 

b5/b4 

Gao, 1996 

Nonnalized Difference 
Wetness Index 

(b2 - b4)/(b2 + b4) 

Crist & Cicone, 1984 

Nonnalized Bum Ratio 

(b4 - b7)/(b4 + b7) 

Key & Benson, 2006 

Tasseled Cap Brightness 

0.3037xbl + 0.2793 x b2 + 0.4743 x b3 + 
0.5585 x b4 + 0.5082 x b5 + 0.1863 x b7 

Johnson & Ross, 2005 

Tasseled Cap Greenness 

-0.2848xbl - 0.2435 x b2 - 0.5436 x b3 + 
0.7243 x b4 + 0.084xb5 - 0.18 x b7 

Johnson & Ross, 2005 

Tasseled Cap Wetness 

0.1509 x b 1 + 0.1973xb2 + 0.3279 x b3 + 
0.3406xb4 - 0.7112 x b5 - 0.4572 x b7 

Johnson & Ross, 2005 


Note: band numbers in the formulas refer to the Landsat 5 band order. 


3.5 Accuracy assessment and the sample size effect 

We conducted both qualitative and quantitative accuracy assessments. First, a confusion 
matrix was created for each individual year, and associated accuracy measures were 
generated, including the overall accuracy, kappa, producer’s accuracy and user’s 
accuracy. Second, we also visually evaluated the results by comparing the mapped spatio- 
temporal pattern of mortality with those observed by the Forest Health Monitoring Aerial 
Detection Survey (ADS) data (Friedl & Brodley, 1997). The ADS data are generated by 
aerial sketchmapping, a technique for observing and mapping forest damage from an 
aircraft. Although ADS had been criticized for its imprecision at fine scales due to 
reasons, such as misregistration and the subjectiveness of the operators (Meigs et al., 
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2011), it is reliable when used at coarse scales (Friedl & Brodley, 1997). A pilot study 
initiated in the Rocky Mountain Region demonstrated that all three types of accuracies 
for the MPB detection in lodgepole pines exceeded 80% when the spatial tolerance was 
increased to 500 m (Friedl & Brodley, 1997). We thus had confidence that the epicenters 
and patterns of spread derived from ADS polygons could be used as a reference dataset to 
compare against our results. We also altered the size of the training sample dataset from 
106 to 1060 in increments of 10% and tested the effects of sample size on accuracy. This 
produced 10 corresponding sets of results to test the effect of sample size on accuracy 
metrics for both the MLC and RF classifications. 

4 Results 

4.1 Parameter calibration for temporal segmentation and decision tree labeling 

Parameter calibration in the forest growth trend analysis involved two parts. At the 
temporal segmentation stage, the set of parameters producing the highest score from all 
144 combinations was chosen (Table 3). By using the same set of training samples, we 
also determined the four threshold values in the DT. We set the allowed variability in 
NBR for stable segments (thre_magl) to be 20 (the unit for all thresholds was NBR x 
100) as reflected from Figure 4. The cut-off NBR value of vertices (thre vertexl) to 
distinguish healthy forest from disturbances was 350. Threshold Thre_mag2, used to 
separate clearcuts from MPB mortality, was set at -150. A higher degree of overlap was 
observed between the post-clearcut and MPB mortality NBR values, and thus, we set 
thre_vertex2 to the third quantile of NBR values in clearcuts, which was 50. 


Magnitude threshold NBR vertex vaule 




Figure 4. Tukey boxplot for (a) the averaged annual normalized bum ratio (NBR, 
multiplied by 100) changes of stable, regrowth, MPB mortality and clearcut plots; 

(b) current vertex values of post-MPB mortality, post-clearcut events and healthy plots. 
The bottom and top of the box are the first and third quartiles, and the band inside the box 
is the median. The whiskers represent the 1.5 interquartile range of the lower and 
upper quartile. 
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4.2 The performance of forest growth trend analysis 

The forest growth trend analysis achieved an overall annual accuracy ranging from 
86.74% to 94.00%, and the average was 90.31% (Table 6). Except for the user’s accuracy 
for MPB mortality (84.81%) and the producer’s accuracy for clearcuts (77.30%), the 
mean value for all types of accuracies was above 85%. Larger standard deviations of 
accuracies among different years were found in the producer’s accuracy for healthy forest 
and the producer’s accuracy for clearcuts. As indicated from the confusion matrix (Table 
S2) and Figure 5, healthy pixels were well distinguished from the two disturbance types 
in the first few years and became less separable from MPB mortality after 2005. Pixels 
with MPB mortality were confused with either healthy or clearcut classes, but the 
proportion was not high. Confusion between clearcuts and healthy classes occurred in the 
first few years of our study and then between clearcuts and MPB mortality in later years. 
Visual examination showed that the shifting of MPB outbreak epicenters and the 
expansion of the outbreak boundaries was congruent between the two independent 
datasets (ADS and validation points). Maps of the onset year of MPB mortality revealed 
that the MPB outbreak in Grand County was initiated at multiple locations and spread 
outward from those locations over time to affect nearly all habitats with suitable hosts 
(Figure 6). 



Figure 5. The proportion of validation samples that were erroneously labeled in other 
classes. Notes: H, healthy; M, MPB mortality; C, clearcut; H2M is explained as “healthy 
samples classified to MPB mortality.” 

4.3 Comparison with single-date classification and sample size effects 

The accuracy achieved by the single-date classification ranged between 37.20% and 
75.20% and was lower than the time series classification, running on either the raw or 
combined spectral bands trained by the same set of samples (Table 7). We also observed 
that the time-series classification accuracies were relatively stable over time. In contrast, 
the single-date image classifications had high year-to-year variability in accuracy, which 
was shown by a distinct, increasing trend of the producer’s accuracy for MPB mortality 
from 2005 to 2011. The producer’s accuracy for healthy forest decreased substantially 
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over the same time period. Table 7 also listed the accuracies generated from RF and MLC 
with only the raw spectral bands and the combined use of raw spectral bands and spectral 
indices. The accuracy of the RF classification improved when the spectral indices were 
included as predictors, and this effect was not evident with the MLC. Between the two 
single-date classification methods, RF had better overall performance than MLC. 
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Figure 6. (a) Classification results of the forest growth trend analysis in the years 2000, 
2005 and 2011; (b) maps of the onset year, duration and magnitude of mountain pine 
beetle (MPB) mortality. 


By varying the number of training samples used to train the MLC and RF classifiers in 
10% increments to a maximum of 1060 training samples, we obtained 10 sets of results 
for the single-date classifications. Figure 7 shows the relationships between accuracy and 
sample size, for overall accuracy generated from the MLC and RF classifiers, as well as 
the producer’s accuracy and the user’s accuracy for the MPB mortality class. Despite 
some small fluctuations, all accuracies increased with the training sample size, and the 
changes were especially pronounced for the initial increases in sample size. The highest 
overall accuracies were achieved when 100%, 100% and 90% of the training samples 
were used by RF (89.43%, 90.48%, 88.66%) and 100% of the training samples were used 
by MLC (88.41%, 82.93%, 85.90%). This improvement was substantial compared to the 
results based on the original set of 106 training samples. The overall accuracy for RF 
increased by 30% on average from the smallest to the highest training sample size, 8% 
for the producer’s accuracy and 33% for the user’s accuracy. For MLC, the average 
increase was 40% in overall accuracy, 31% in the producer’s accuracy and 43% in the 
user’s accuracy between the smallest and largest training sample sizes. We also noticed 
that after a certain training sample size, the accuracy stopped increasing or even slightly 
decreased, and the saturation position was around the 40% sample size (or 424 points) for 
all accuracy metrics and both classifiers. Accuracies after the saturation point were 
comparable to the time series classification results. 

Table 6. Accuracy assessment results for individual years of the time series 
classification. 


Year 



Healthy 

MPB Mortality 

Clearcut 


OA 

Kappa 

PA 

UA 

PA 

UA 

PA 

UA 

2000 

89.00 

82.21 

100.00 

83.02 

88.89 

94.12 

65.00 

100.00 

2001 

90.00 

84.23 

100.00 

95.74 

100.00 

79.49 

58.33 

100.00 

2002 

89.00 

82.88 

100.00 

95.56 

93.75 

78.95 

64.00 

94.12 

2003 

93.00 

88.92 

100.00 

93.18 

94.87 

92.50 

75.00 

93.75 

2004 

92.00 

87.41 

100.00 

100.00 

89.74 

89.74 

78.95 

78.95 

2005 

91.00 

85.51 

91.49 

97.73 

97.14 

80.95 

77.78 

100.00 

2006 

86.74 

78.99 

72.97 

96.43 

95.24 

78.43 

94.74 

94.74 

2007 

94.00 

90.49 

92.11 

100.00 

97.67 

89.36 

89.47 

94.44 

2008 

91.00 

86.20 

82.35 

96.55 

97.50 

84.78 

92.31 

96.00 

2009 

91.00 

85.81 

87.10 

96.43 

100.00 

84.62 

80.00 

100.00 

2010 

90.00 

84.56 

97.14 

97.14 

95.00 

84.44 

72.00 

90.00 

2011 

87.00 

79.64 

81.25 

92.86 

95.35 

80.39 

80.00 

95.24 

mean 

90.31 

84.74 

92.03 

95.39 

95.43 

84.81 

77.30 

94.77 

SD 

2.18 

3.45 

9.25 

4.46 

3.47 

5.47 

11.35 

5.92 


Notes: OA, overall accuracy; PA, producer’s accuracy; UA, user’s accuracy; SD, 
standard deviation. 
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Table 7. Accuracy comparison among the time series classification, random forests 
(RF) and maximum likelihood classifier (MLC) with the same set of 2005, 2009 
and 2011 validation samples. Only the mountain pine beetle (MPB) mortality class 
is included here. 


Year 


Time-Series 
Trend Analysis 


RF 

MLC 


UA 

80.95 

raw * 
68.00 

combined 

67.30 

raw 

71.96 

combined 

88.00 

2005 

PA 

97.14 

18.78 

62.40 

85.08 

65.90 


OA 

91.00 

63.82 

74.40 

75.20 

73.50 


UA 

84.62 

46.15 

100.00 

39.63 

0.00 

2009 

PA 

100.00 

35.43 

48.70 

100.00 

100.00 


OA 

91.00 

53.47 

72.90 

40.94 

38.70 


UA 

80.39 

78.47 

77.80 

37.92 

22.20 

2011 

PA 

95.35 

83.06 

10.80 

100.00 

100.00 


OA 

87.00 

66.62 

40.40 

38.15 

37.20 


Notes: “raw” denotes classification with solely spectral data, and “combined” denotes 
classification with the combined use of spectral and index data; PA, producer’s 
accuracy; UA, user’s accuracy; OA, overall accuracy. 



0 20 40 60 80 100 0 20 40 60 80 100 0 20 40 60 80 100 


OA PA UA 

Figure 7. Relationships between accuracy and sample size for the maximum likelihood 
classifier (MLC) and random forests (RF) using multiple training sets. Graphs arranged 
from left to right display the trends of overall accuracy, the producer’s accuracy and the 
user’s accuracy for the mountain pine beetle (MPB) mortality class as the training sample 
proportion increases from 0.1 to one. Blue and solid lines represent RF results, whereas 
green and dashed lines represent MLC results. Red arrows indicate the mean accuracies 
produced by the forest growth trend analysis in the years 2005, 2009 and 2011. 


5 Discussions 
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For the single-date classification, we found the following: (1) The RF achieved better 
overall accuracy than MLC in detecting forest disturbances. MLC relies on assumptions 
about the data distribution (e.g., normally distributed), whereas the ensemble learning 
techniques used by RF do not (Ghimire et ah, 2010). Meanwhile, by grouping many 
individual classifiers, RF combined the strengths and minimized the weaknesses of each 
(Pal & Mather, 2003). (2) MLC is less affected by the selection of input layers than RF. 
RF relies on all outputs produced by each individual tree that is built upon the random 
selection of input layers, and thus, a larger amount of features could create a higher 
probability of “irrelevant” features being selected; but, when this occurs, it is usually 
negated by using the ensemble of results from individual trees. (3) The accuracies of 
single-date image classifications showed prominent year-to-year variability. The 
omission error for MPB mortality declined consistently from 2005 to 2011, whereas the 
omission error for healthy forest increased significantly. During this period, Grand 
County was in a steady transition from a landscape dominated by healthy forests to a 
landscape dominated by “grey-stage” forests. Not coincidently, the declining trend in the 
proportion of healthy forest matches the declining tendency of its omission errors over 
the years. 

From the experiment on the sample size effect, on the one hand, we could detennine the 
optimal sample size in the single-date MPB mortality mapping, which, to our knowledge, 
has not been discussed before. Both MLC and RF were sensitive to sample size, and 
several studies are in accordance with our conclusion (Foody, 2009; Li et ah, 2013; Pal, 
2005). Parametric algorithms are well known for their requirement of representative 
samples, and this requirement is usually satisfied by having a sufficiently large sample size 
(Rodriguez-Galiano et ah, 2012a). Studies on sample-size effects for non-parametric 
classifiers are relatively rare. Rodriguez-Galiano et al. (2012b) observed a significant 
reduction in overall accuracy and the kappa values when the sample size reduced more 
than 50% of a non-parametric RF applied to Landsat images. 

In our case, if we take both accuracy and labor/economic/time costs into consideration, 
the most optimum sample size for single-date classification is four to five times the 
original training sample size, which included 106 samples. On the other hand, we could 
compare the efficiency of single- and multi-date classifications. By using the same 106 
training samples, our time series approach substantially outperfonned the single-date 
classification, and their accuracies were similar only after expanding the training sample 
size to a large number (e.g., 10 times), which matches the finding from another study 
concluding that either single- or multi-date methods in Landsat-based forest disturbance 
mapping could result in high classification accuracy (Meddens et al., 2013), but our 
results confirm their findings only for large training sample sizes. 

Although the forest growth trend analysis method is promising in many ways, several 
shortcomings exist: (1) Regardless of the high overall accuracy, the commission error for 
clearcuts and the omission error for beetle mortality were relatively high. The 
commission errors in clearcuts could be caused by the spectral similarity between post- 
clearcut pixels dominated by successional vegetation and healthy forest pixels. Errors of 
commission for MPB mortality usually happened in areas hosting a mixture of grey stage 
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trees and bare land, which could lead to a misclassification to clearcuts. (2) Parameter 
calibration was the most time-consuming part in the whole procedure and requires the 
expertise of image analysts. In some sense, the identification of various disturbance types 
using temporal trajectories is more challenging than the visual interpretation of false- 
color or true color image composites. Several potential solutions include the use of a tool 
called TimeSync, which was developed to aid with the temporal segmentation calibration 
and to verify the output from Landtrendr (Cohen et al., 2010). In addition, the Forest 
Inventory and Analysis (FIA) Program of the U.S. Forest Service consistently collects 
information on status and trends in forest species, health status, in total tree growth, 
mortality and removals by harvest. If FIA plots, or similar types of data, and their exact 
locations are available to users, they could aid the identification of the plots’ disturbance 
history. In terms of the labeling, although the manmade decision tree, as shown in Figure 
3, can be easily interpreted, a machine learning-based decision tree could be employed to 
make this step fully automatic. 

6 Conclusions 

Reconstructed forest disturbance histories are needed for a range of applications, and our 
paper tested an effective and efficient approach to generating disturbance histories using 
the Landsat time series stack. Our algorithm was designed to take advantage of 
differences in the magnitude and direction of changes in spectral indices, which are 
proxies for the physiological processes following different types of forest disturbances. 
The site-specific accuracy assessment in the three validation years indicated that the 
overall accuracy ranged from 86.74% to 94.00%, and the mean value of all types of 
accuracies was above 85%. Especially, our method outperformed standard classification 
approaches (maximum likelihood classification and random forest) by 15.8% to 52.3% 
when small training sample sizes (106 samples per study area) were used, which is a cost¬ 
saving advantage. Another strength of our strategy is that most parameters, especially at 
the attributing stage, are easily interpreted and modified, which will especially benefit 
users tackling similar problems of labeling change detection results. 

In conclusion, our results demonstrate the potential for a forest growth trend analysis to 
characterize forest disturbance events with an optimal sampling scheme, both in 
economic and time terms, at a satisfying mapping accuracy level. Some short-term goals 
for our further studies will include improving the labeling step by replacing it with an 
automated machine learning technique, linking the magnitude of change with more 
ecologically meaningful properties, such as vegetation cover, extending the algorithm to 
the ecoregion scale, building spatio-temporal dynamic outbreak models on the time series 
maps with projections into the future. In addition to those applications, the spatially 
explicit infonnation produced by our approach could also serve as raw material to 
facilitate other research to uncover new patterns of disturbances occurring on the 
landscape and to investigate the drivers generating the patterns of disturbances. 
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Table SI. MPB host species in LF Refresh 2001. Their vegetation type code, vegetation type, system group name, canopy cover, area 
of occupancy in Grand County, associated omission error and commission error are listed. 


Code 

Vegetation type 

System group name 

Canopy 

cover 

Area 

(km 2 ) 

Omission 
error (%) 

Commission 

error(%) 

2049 

Rocky Mountain Foothill Limber Pine- 
Juniper Woodland 

Limber Pine Woodland 

Open 

0.73 

NA 

NA 

2050 

Rocky Mountain Lodgepole Pine Forest 

Lodgepole Pine Forest and Woodland 

Closed 

140.01 

54.50 

18.60 

2051 

Southern Rocky Mountain Dry- 
MesicMontane Mixed Conifer Forest 
and Woodland 

Douglas-fir-Ponderosa Pine- 
Lodgepole Pine Forest and Woodland 

Closed 

43.58 

53.30 

46.20 

2057 

Rocky Mountain Subalpine-Montane 
Limber-Bristlecone Pine Woodland 

Limber Pine Woodland 

Open 

2.25 

NA 

NA 

2117 

Southern Rocky Mountain Ponderosa 
Pine Savanna 

Ponderosa Pine Forest, Woodland 
and Savanna 

Sparse 

0.004 

NA 

NA 


Note: NA - No data in the LANDFIRE national quality assessment report. 


UJ 

UJ 






Table S2. Confusion Matrix of the forest growth trend analysis for each year from 2000 
to 2011. 


Reference data 


2000 

H 

A 

C 

UA(%) 

H 

44 

4 

5 

83.02 

Classification A 

0 

32 

2 

94.12 

C 

0 

0 

13 

100.00 

PA(%) 

100.00 

88.89 

65.00 



2001 

H 

A 

C 

UA(%) 

H 

45 

0 

2 

95.74 

A 

0 

31 

8 

79.49 

C 

0 

0 

14 

100.00 

PA(%) 

100.00 

100.00 

58.33 


2002 

H 

A 

C 

UA(%) 

H 

43 

1 

1 

95.56 

A 

0 

30 

8 

78.95 

C 

0 

1 

16 

94.12 

PA(%) 

100.00 

93.75 

64.00 


2003 

H 

A 

C 

UA(%) 

H 

41 

1 

2 

93.18 

A 

0 

37 

3 

92.50 

C 

0 

1 

15 

93.75 

PA(%) 

100.00 

94.87 

75.00 


2004 

H 

A 

C 

UA(%) 

H 

42 

0 

0 

100.00 

A 

0 

35 

4 

89.74 

C 

0 

4 

15 

78.95 

PA(%) 

100.00 

89.74 

78.95 


2005 

H 

A 

C 

UA(%) 


34 



H 

43 

1 

0 

97.73 

A 

4 

34 

4 

80.95 

C 

0 

0 

14 

100.00 

PA(%) 

91.49 

97.14 

77.78 


2006 

H 

A 

C 

UA(%) 

H 

27 

1 

0 

96.43 

A 

10 

40 

1 

78.43 

C 

0 

1 

18 

94.74 

PA(%) 

72.97 

95.23 

94.74 


2007 

H 

A 

C 

UA(%) 

H 

35 

0 

0 

100.00 

A 

3 

42 

2 

89.36 

C 

0 

1 

17 

94.44 

PA(%) 

92.11 

97.67 

89.47 


2008 

H 

A 

C 

UA(%) 

H 

28 

0 

1 

96.55 

A 

6 

39 

1 

84.78 

C 

0 

1 

24 

96.00 

PA(%) 

82.35 

97.50 

92.31 


2009 

H 

A 

C 

UA(%) 

H 

27 

0 

1 

96.43 

A 

4 

44 

4 

84.62 

C 

0 

0 

20 

100.00 

PA(%) 

87.10 

100.00 

80.00 


2010 

H 

A 

C 

UA(%) 

H 

34 

0 

1 

97.14 

A 

1 

38 

6 

84.44 

C 

0 

2 

18 

90.00 


35 



PA(%) 97.14 95.00 72.00 


2011 

H 

A 

C 

UA(%) 

H 

26 

1 

1 

92.86 

A 

6 

41 

4 

80.39 

C 

0 

1 

20 

95.24 

PA(%) 

81.25 

95.35 

80.00 
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Chapter 3 Forest disturbance regimes, interactions, and successional 
pathways in a mountain pine beetle infected ecoregion 

Abstract 

The pine forests in the southern portion of the Rocky Mountains are a heterogeneous 
mosaic with some areas receiving extensive and intensive stress from human activity, fire, 
and mountain pine beetles (MPB; Dendroctonus ponderosae ), and other areas in various 
stages of recovery from disturbances. Understanding the underlying disturbance regimes, 
their interactions, and disturbance-succession pathways are crucial for adapting 
management strategies to mitigate their impacts. With this purpose, we depicted the 
forest disturbance and recovery history over a 13-year time series in the Southern Rocky 
Mountains Ecoregion with Landsat image stacks, aided by a novel sampling strategy that 
takes advantage of the scene-overlap among adjacent Landsat images, and an automated 
workflow integrating a temporal segmentation technique developed from Landsat-based 
Detection of Trends in Disturbance and Recovery algorithm, a set of spatio-temporal 
segment features and Random Forest classifier. The mean overall accuracy for all scenes 
was 82%. Analysis of the disturbance trends revealed that MPB contributed to 70% of the 
total disturbed area. We found that burn severity was largely unrelated to the severity of 
pre-fire beetle outbreaks in this region, where the severity of post-fire beetle outbreaks 
has been shown to decrease with bum severity. Approximately half the clearcut and 
burned areas were in various stages of recovery, but the regeneration rate was much 
slower for MPB-disturbed sites. Pre-fire beetle outbreaks and subsequent fire produced 
positive compound effects on seedling re-establishment in SRME, which suggests that 
serotinous lodgepole pine is resilient to wildfires occurring after beetle mortality. 

1 Introduction 

The vegetation of the Southern Rocky Mountains Ecoregion (SRME) is a heterogeneous 
disturbance-derived mosaic, where fire, blow-down, insect outbreaks and clearcut are 
regarded as the primary disturbance types (Barbour & Billings, 2000; Veblen et al., 1994). 
Under the pressure from warming climate and human interventions, the disturbance and 
recovery regimes, which refer to the spatial and temporal DR dynamics over a long time 
period (Turner, 2010), are undergoing rapid changes in the subalpine zones of Rocky 
Mountains. Fire is perhaps the most commonly studied disturbance agent (eg. Justice et 
al., 2002; Simon et al., 2004; Stroppiana et al., 2000) and the total area burned and the 
number of large fires have risen dramatically across much of the western United States 
over the past 25 years (Calkin et al., 2005; Stephens, 2005). Meanwhile, massive conifer 
mortality since the mid-1990s in this ecoregion has been related to outbreaks of bark 
beetles: primarily mountain pine beetle (MPB), spruce beetles, and western balsam bark 
(Smith et al., 2015). Among them, MPB mortality has intensively expanded in high 
elevation conifers in the Rocky Mountains reaching epidemic levels not previously 
recorded (Gibson et al., 2008). The amount of carbon in trees killed by fires and by MPB 
outbreaks has been quite similar over the past decade (Hicke et al., 2013). 
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The subsequent recovery of forest vegetation is also highly spatially and temporally 
variable, since disturbances couple with environmental factors to determine the rates and 
pathways of succession (Pfulgmacher et ah, 2013). Both disturbance and recovery (DR) 
fundamentally shape the trajectory of forest carbon storage, influence forest structure, 
composition, biomass, land surface water and radiation dynamics, and affect feedbacks 
between forest ecosystems and climate (Baldocchi, 2008; Bonan, 2008; Franklin et ah, 
2002; Gough et ah, 2007; Kurz et ah, 2008). Thus, understanding DR patterns, as well as 
their interactions are critical for assessing their impacts on the forest resilience 
mechanisms, and benefit relevant studies, such as habitat forecasting, general circulation 
modeling, international carbon pricing negotiation, and biodiversity conservation (Gibbs 
et ah, 2007; Liang et ah, 2014a; Running, 2008). Under this context, the first goal of this 
paper is to update the knowledge on the historical and current patterns of DR regimes in 
SRME, with our focus on the changes of MPB disturbance over the past decade. 

The occurrences of multiple disturbances in this ecoregion provide us a good opportunity 
to investigate their interactions - a raising concern in ecology because they may cause 
ecosystems to be pushed beyond thresholds (White & Pickett, 1985; Groffman et ah, 

2006; Turner, 2010), yet a challenging question due to the scarcity of quantitative data. A 
disturbance may be linked with another by changing its extent, severity, or probability of 
occurrence. For SRME, there is substantial interest in how the extensive outbreaks of 
bark beetles may affect future wildfire (Turner, 2010). While spruce beetle outbreaks 
were detected to have no overlap with stand-replacing fires in the Colorado Rocky 
Mountains (Veblen et ah, 1994), conclusions regarding MPB and fire interactions are 
more controversial. Increased (Hopkins, 1909; Schoennagel et ah, 2012), decreased 
(Simard et al., 2011), and unrelated (Bond et ah, 2009; Veblen et al., 1994) active crown 
fire potentials were found following MPB outbreaks, depending on the outbreak stage 
and burning conditions (Harvey et al., 2014). Additionally, there is no consensus whether 
MPB outbreaks following fires affect forest vegetation. Although it is presumed that fire- 
injured forests are more susceptible to beetles (Geiszler et al., 1984), some have pointed 
out that wildfire is an unlikely driver of outbreaks by MPB (Powell et al., 2012). Till now, 
no study has been conducted to investigate their interactions across the whole South 
Rockies. Accordingly, the second goal we attempted to address was to determine whether 
MPB outbreaks and wildfire are linked disturbances in this ecoregion. 

If such linkages exist between disturbances, they may produce compound effects on 
ecosystems, which could be reflected in amplified severity beyond that of any singular 
disturbance, or the creation of novel characteristics that would not normally happen in 
isolation (Buma & Wessman, 2011), such as the shifting of dominant forest cover types 
(Johnstone et al., 2010), and reduced regeneration rates (Buma & Wessman, 2011). 
Understanding the consequences of those interactions is critical for conservation and 
management purposes. Our third goal is thus focused on revealing the effects of various 
compounded disturbances on lodgepole pine forest regeneration in SRME as compared 
with singular disturbances. 

Long-term data of historical and current disturbances that has same spatial and temporal 
extent are lacking, but required to address our 3 research goals. Dendrochronology is the 
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most applicable and prevalent technique, but can be limited in geographical scope and 
may not detect all disturbance types. Forest Health Monitoring Aerial Detection Surveys 
is the major state and federal forest health monitoring effort since the mid-twentieth 
century (Man, 2010), but it has been criticized for spatial imprecision caused by 
misregistration as well as operators’ subjectiveness, and the temporal shift between the 
attack onset and the spectral signature changes of canopy (Meigs et al., 2011; Kautz, 
2014). Satellite spectral indices have been proven to be useful in accurately delineating 
forest disturbances over long time spans at broad spatial scales (Wulder et ah, 2006), but 
if done with single date imagery, relative subtle changes arising from insect attacks are 
more difficult to monitor than prominent changes associated with stand replacing events 
(Skole & Tucker, 1993). We recently completed a twelve-year forest disturbance 
mapping in Grand County, Colorado with good success using Landsat time-series stack 
(Liang et ah, 2014b). But so far, this approach has not been applied over large scales due 
to its dependency on human interpretation in attributing disturbance events. Here, we will 
advance this approach by developing an efficient and effective sampling strategy and 
automating the event attribution step, to serve the purpose of decadal forest DR mapping 
at the ecoregion scale. 

2 Study area 

Our study area in the SRME extends south from southern Wyoming through Colorado to 
the northern New Mexico (lower left: ~35°N, 109°W; upper right: ~43°N, 103°W) 

(Figure 1). To avoid confounding different types of insect mortality (e.g. spruce vs. pine 
beetles), we confined our analysis to the extent of MPB host species in the LANDFIRE 
Refresh 2001 Existing Vegetation Type data layer that was developed with circa 2001 
Landsat imagery (Nelson et al., 2013), and was smoothed with a 3 X 3 majority filter to 
reduce noise and minimize the reported high omission errors (National c2001 Assessment; 
Figure 1). Pine forest selected as MPB host species included lodgepole pine (P. contorta; 
95.4%), limber pine ( P.flexilis ; 2.9%) and ponderosa pine ( P. ponderosae; 1.7%). Their 
total coverage was approximately 60,675 km 2 . Hereafter, these forest types are referred to 
as ‘forest’. 
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Figure 1. The coverage of Southern Rocky Mountains Ecoregion (U.S. Environmental 
Protection Agency, 2013), and the frame of 14 intersected Landsat scenes with trimmed 
edges. 
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3 Methods 


3.1 Remote sensing image pre-processing 

Landsat time series stacks (LTSS) are well suited to this study, because of their free 
availability, long time period coverage, broad spatial extent, multispectral characteristics, 
and temporal continuity. All Landsat images were acquired from the Landsat Surface 
Reflectance Climate Data Record generated by the Landsat Ecosystem Disturbance 
Adaptive Processing System (URL: http://earthexplorer.usgs.gov; accessed: 22 August 
2014). We restricted our analysis to images acquired June through October, during which, 
June to September is the most active growing season for vegetation and thus provides the 
strongest spectral contrast between the healthy and beetle mortality forest. Images in 
October usually have better quality than summer in this mountainous region and were 
used as substitutions when clouds heavily contaminated all images in the ideal time 
window. We tried to find images containing less than 20% cloud cover, and if no such 
data existed, two or more cloudy images with similar phenology characteristics were used 
as replacements to guarantee one clear observation of a MPB infected area for most 
regions. Fourteen Landsat scenes intersected the boundary of SRME, and 247 images 
were collected from 1999 to 2011 in total with approximately 18 images per scene 
(Figure 2). We used the Normalized Bum Ratio (NBR; Key & Benson, 2005) to 
reconstruct the spectral trajectory, because of its sensitivity to disturbance and its capacity 
to minimize radiometric differences among images (Kennedy et al., 2010). 



Figure 2. Scene acquisition dates of the fourteen Landsat scenes by year. Supplementary 
images used to fill the cloud gaps of the image in the same year were marked as magenta 
stars. 


41 







3.2 Reference sample selection 

To improve the sampling efficiency at the ecoregion scale, we developed a new strategy 
by intensively selecting samples in the Landsat scene-overlapping area (SOA) between 
images along the same row, which occupies 20% to 40% coverage of the entire scene. 

The intention is to save efforts in image interpretation or on-site data collection. 
Meanwhile, although 20%-40% of the full scene coverage is a good representation of the 
landscape, to make the final sample set an unbiased indication of what the population is 
like, we then automatically add spectrally dissimilar samples from the non-overlapping 
area (NOA). Specifically, our sampling strategy involved the following steps (Figure 3): 

1) Clustering. The LTSS was divided into SOA and NOA, and K-Means clustering 
was applied on them separately, in order to group pixels sharing similar 
trajectories together. 

2) Calculate spectral distances between cluster centers. The spectral distance 
between the center of each SOA and NOA cluster was calculated. 

3) Remove spectrally similar NOA clusters. NOA clusters that shared similar spectral 
signature with any SOA cluster would be treated as duplicated samples and 
removed. That is, only NOA clusters with their centers to all SOA cluster centers 
larger than one distance standard deviation of the distance set are retained. The 
combination of all SOA clusters and non-duplicated NOA clusters was thus able 
to represent the majority spectral space of the overlapped images. 

4) Stratified random sampling and event labeling. Finally, we conducted stratified 
ransom sampling on the combined cluster set. For each sample, we interpreted 
their DR history with the aid of NBR trajectories, false-color composite Landsat 
layers, and high resolution Google Earth images. 
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Figure 3. Sampling work flow. SOA: Landsat scene-overlapping area; NOA: non¬ 
overlapping area; dis - spectral distance; Cl: one distance standard deviation of the 
spectral distance set. 


3.3 Temporal segmentation 

Temporal segmentation is one approach used to extract useful information from spectral 
trajectories, and has been successfully implemented with Landsat time-series in Landsat- 
based detection of Trends in Disturbance and Recovery tool (Landtrendr; Kennedy et ah, 
2010). The idea is to decompose the often-noisy time series into a sequence of straight- 
line segments that can represent distinct vegetation states. The general procedure is to 
first remove noise-induced spikes, then identify potential vertices that can separate 


43 












































distinct periods before and after them, apply a set of trajectory-fitting algorithms to best 
capture the salient events happening in the trajectory, and finally choose the best model 
using an F-statistic (Kennedy et ah, 2010). The quality of temporal segmentation is 
controlled by a group of parameters, and the calibration method was outlined in Liang et 
al. (2014a). For each pixel, the outputs include several fitted segments with a set of 
features describing the trajectory. 

3.4 Segment-based random forest classification 

We used temporal segments as the basic unit in the Random Forest (RF) classification 
instead of pixels. Because different DR events have distinct segment features, such as 
duration, severity, and onset, the combination of which can thus result in better 
differentiation. For instance, clearcuts could be recognized as segments with rapid and 
severe spectral value changes, whereas the changes of MPB mortality segments tend to 
be long, but slow. Besides the features of the focal segment, we also constructed a set of 
spatial neighborhood features and temporal neighborhood features to feed the RF 
classifier, which were designed based on the assumption that events that happened within 
a proximal zone at the same time or to a fixed location at different times are highly 
correlated with each other. The technical details of feature construction are explained as 
following: 

Segment features : Vertex year is the starting year of the event occurrence. The fitted 
vertex value is the NBR value that comes out of Landtrendr segmentation. Magnitude is 
defined as the absolute NBR change between the head and tail vertex of the focal 
segment (ANBR — NBR head — NBR taii ). Duration denotes the time span from the head 
to the tail vertex, and slope is the ratio of magnitude to duration (Figure 4). 



Figure 4. An illustration of the segment features and the neighborhood features listed in 
Table 1. 
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Temporal neighborhood features included the duration, magnitude, and slope of the 
prior and posterior segments of the focal one. 

Spatial neighborhood features were calculated based on the slopes and vertex values of 
the eight spatially adjacent pixels over the duration of the focal segment. The slope set is 
denoted as N s and the vertex value set is defined as Nf: 

Ns {(SaljSblj • • • 5 Snl), (Sa2,Sb2, • • • ,Sn2), • • •, (Sa8,Sb8, • • • ,Sn8) 

Nf= {(Fal,Fbl,...,Fnl), (Fa2,Fb2,...,Fn2),..., (F a 8,Fb8,... ,F n 8) 

Where, S is slope and F is the fitted vertex value. The first subscript a to n indicates the 
first year a of the time span to the last year n, and the second subscript 1 to 8 represents 
the first to the eighth neighbor pixel. For example, S a i is the slope of the first neighbor 
pixel in the first year of the focal segment duration. 

Following that, we calculated the maximum value for each subset in N s , to highlight the 
event with the most significant change (max8 slope), and computed the mean and standard 
deviation of max8slope to represent the overall condition in the neighborhood system 
(avg_max8slope; std_max8slope). The same processing routine was applied to the Nf, 
with an additional minimum value calculation. 

In total, forty-two segment features falling into fifteen categories were created. Five 
broad classes that represent major forest DR status in SRME were classified: healthy, 
MPB mortality, fire, clearcut, and regrowth. For sudden canopy removal events, we 
separated fires from clearcuts using the U.S. Geological Survey Monitoring Trends in 
Bum Severity database (Eidenshink et ah, 2007). We then applied a plurality vote to the 
ten RF outcomes generated from different training partitions, and produced the annual 
time-series DR maps by applying the segment label to their corresponding years. Finally, 
the classification performance was evaluated via 10-fold cross validation. The overall 
accuracy (OA) was calculated by summarizing the percentage of segments that were 
labeled the same by the RF classifier and our samples, with the assumption that the 
temporal segmentation step is flawless. This assumption is made because it is difficult to 
obtain historical, temporally-explicit ground truth data covering the whole ecoregion. 

3.5 Disturbance regime, interactions and disturbance-recovery pathways 

We derived size, onset year, magnitude, and duration as the key disturbance regime 
parameters for MPB mortality, fires and clearcuts. Onset represents when events are 
initiated over time. Magnitude is an indicator of the percent vegetative cover change 
(Kennedy et ah, 2010), and duration suggests the most dominant event over the last 
decade for a particular location. 

Two types of interactions exist between MPB outbreaks and fire in terms of their 
occurring order: MPB+fire (fire happens after beetle outbreaks), and fire+MPB (beetle 
outbreaks happen following fire). We attempt to detennine their linkages by investigating 
how recent bark beetle outbreaks influenced subsequent bum severity, and whether burn 
severity affects the severity of pre-fire beetle outbreaks. To answer the question, we first 
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allocate pixels into four groups based on their disturbance history. MPB+fire and 
fire+MPB target groups refer to the assemblage of pixels containing the process of the 
two interaction modes. MPB / fire population group is the collection of pixels not 
showing the interaction patterns but possessing beetle outbreaks / fire events. The 
inclusion of the population groups was designed to control for the interaction groups to 
be compared with a spectrum of outbreak or bum severity occurring over the study region. 

We then assessed the magnitude values for MPB events of the two target groups against 
those of the MPB population group separately using analysis of varianace (ANOVA), and 
the same procedure was applied to fire. The significance test on the means of the target 
and population groups provided a way of examining the effects of one disturbance event 
on the severity of the other. For instance, if no significant difference was found between 
the MPB outbreak magnitude in the fire+MPB target group and the MPB population 
group, the fire probably did not play a role in changing the severity of beetle outbreaks. 

Finally, we depicted the disturbance-recovery pathway by summarizing the area size and 
proportion of each disturbance event that transitioned to either another disturbance type 
or to the recovery stage. 

4 Results 

4.1 Sampling and mapping 

To demonstrate how the new sampling strategy functions, we displayed the centers of 
SOA and NOA clusters from one scene as an example in the Tasseled Cap space for 
visualization purpose, so that the concentrations of data within the total data volume can 
be most readily observed (Figure 5). It was intuitively evident that the centers of the 
duplicated NOA clusters were highly overlapped with the SOA clusters, whereas those of 
the few non-duplicated NOA clusters were much further away from them. This 
demonstrated that our sampling strategy could promote both the efficiency and accuracy 
in choosing representative samples when compared with random sampling over the entire 
scene. 

The inclusion of spatio-temporal features has been proven to be effective as revealed 
from our accuracy evaluation results. The segment based OA for individual scenes varied 
from 74% to 92.5%, with the mean value at 85.9%, and if no such input is considered, the 
classification will yield unifonnly lower accuracies for all scenes, with the mean value at 
82%. 
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Figure 5. Centers of SOA, NOA and duplicated NOA clusters in the Tasseled Cap space 
of one Landsat scene (path: 32; row: 34). Each axis represented the absolute changes in 
brightness, greenness, and wetness over the study period. 

4.2 Disturbance regime descriptors 

Spatial distribution dynamics: Mountain pine beetle outbreaks had a clustered pattern 
(Figure 6). Early in the time span of this study, MPB outbreaks were scattered across the 
SRME and lasted for a limited number of years. However, by the early 2000’s, the MPB 
epidemic in northern Colorado and southern Wyoming along the continental divide 
appeared and persisted over several years (Figure 6). 
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Figure 6. Spatial pattern of onset year (a), NBR change magnitude (b) and outbreak 
duration (c) for MPB disturbance. 
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Disturbance size: During the study period, fires, MPB mortality, and clearcuts together 
affected 10,504 km 2 of forest in this region (17.3% of the host forest area). Mountain pine 
beetle mortality contributed most to the forest loss (7,359 km 2 ). The disturbance event 
affecting the second largest area was fire (2,467 km 2 ). Clearcut events accounted for 678 
km 2 of disturbed area. There were 2,539 km 2 area classified as forest recovering from 
disturbances (Figure 7). 



Figure 7. Area of forests in mountain pine beetle (MPB) mortality, clearcut, fire and 
regrowth classes from 1999 to 2011 in Southern Rocky Mountain Ecoregion. 


Onset time: Time series of the disturbed area revealed a constantly increasing trend of 
MPB caused tree mortality, and two apparent peak years were observed, separately 2001 
and 2008 (Figure 8a). From 2005, the amount of area with new MPB mortality 
accelerated. Notable years with large burned areas (Figure 7) included 2000, 2002 and 
2011, and the other years had a fairly small amount of burned forest area. Clearcut events 
(Figure 7) had a relative stable trend over the time span of this study. The amount of area 
classified as regrowth in the SRME was low in the first five years of the study, but after a 
sharp increase in 2003, the annual area in regrowth stabilized, but was almost twice as 
large as it was in the early part of our study. 

Magnitude: We divided the disturbance magnitude into four levels according to the NBR 
changes, and defined them as low (ANBR:<100), low-middle (ANBR:100 to 200), middle 
(ANBR:200 to 300), and high severity (ANBR: >300). The magnitude of MPB outbreaks 
tended to increase from south to north (Figure 4). More than 30% of the MPB disturbed 
regions were characterized by low-severity beetle mortality (Figure 8b), which 
corresponds to 10.6% canopy cover change per site as inferred from the delta model in 
Kennedy et al. (2010). The most severe outbreaks were observed to happen in northern 
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Colorado and Southern Wyoming, where beetles’ presence became more prevalent in the 
later years (Figure 6). 

Duration: About 20% of beetle outbreaks in this region lasted for 3 years, which 
matched well with the time frame of the commonly acknowledged “green-red-gray” 
stages (Figure 8c). The spatial visualization of duration pattern also revealed that while 
the epicenters of MPB mortality harbored longer-lasting outbreaks (>3 years), their 
peripheral regions were characterized by shorter duration attacks. 


ATree cover(%) 



2000 2002 2004 2006 2008 2010 2012 <100 150 200 250 300 350 400 450 500 >500 0 1 2 3 4 5 6 7 8 9101112 


onset year ANBR(*1000) duration(y) 

Figure 8. The percentage of mountain pine beetle mortality pixels for each (a) onset year, 
(b) nonnalized bum ratio (NBR) / tree cover percentage change interval, and (c) duration. 
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Table 1. Segment features and their definitions. 


Category 

Vertexl_year 

VertexlNBR 

Vertex2_NBR 

Dur 

Mag 

Slope 

max8slope 

avg_max8slope 

std_max8slope 

max8fit 

avgniaxSfit 

std_max8fit 

min8fit 

av g_niin8fit 
std min8fit 


#features Definition 

1 Year of the head vertex 

1 NBR value of the head vertex 
1 NBR value of the tail vertex 

3 Duration of the segment: {precur, curdur, post dur} 

3 Magnitude of change in NBR between the head and tail vertex: {pre mag, cur mag, post mag} 

3 Magnitude/duration: {pre slope, cur slope, post slope} 

8 Max slopes over the duration of the focal segment for its eight spatially adjacent pixels: {max(S a i,Sbi,... ,S n i), 
max(Sa 2 ,Sb 2 ,.. . ? Sn 2 ),.. max(S a 8,Sb8 ? •. *?Sn8)} 

1 Mean ofMax8slope: mean {max(S a i,Sbi,...,S n i), max(S a 2 ,Sb 2 ,.--,Sn 2 ),..., max(S a 8,Sb8,...,S n 8)} 

1 Standard deviation of max8slope: Std {max(S a i,Sbi,.. .,S n i), max(Sa 2 ,Sb 2 ,.. .,Sn 2 ),..max(S a 8,Sb8,.. .,S n s)} 

8 Maximum fitted NBR value over the duration of the focal segment for its eight spatially adjacent pixels: 

{max(F al ,F b i,.. .,F„i), max(F a2 ,F b 2 ,.. - ,Fn 2 ),..max(F a8 ,F b 8,.. .,F„8)} 

1 Mean of max8fit: avg{max(F a i,F b i,...,F n i), max(F a 2 ,F b 2 ,...,Fn 2 ),.--, max(F a 8,F b 8,...,F n s)} 

1 Standard deviation of max8fit: Std {max(F a i,Fbi,...,F„i), max(F a 2 ,Fb 2 ,...,Fn 2 ),..max(F a8 ,F b 8,...,Fn8)} 

8 Minimum fitted NBR value over the duration of the focal segment for its eight spatially adjacent pixels: 

{min(F a i,F b i,.. .,F n i), min(F a2 ,F b 2 ,.. - ,Fn 2 ),..min(F a8 ,F b 8,.. .,Fns)} 

1 Mean of min8fit: avg{min(F a i,F b i,...,F n i), min(F a2 ,F b 2 ,...,Fn 2 ),.--, min(F a 8,F b 8,...,Fn8)} 

1 Standard deviation of min8fit: Std {min(F a i,Fbi,...,F n i), min(F a 2 ,Fb 2 ,..-,Fn 2 ),..., min(F a s,Fb8,...,F n s)} 


Note: all the notations used in sets follow the definition in figure 4. 






4.3 Disturbance interactions and successionalpathways 

Despite the fact that the SRME experienced a diverse range of disturbances, we did not 
observe a significant amount of secondary disturbance happening after a primary 
disturbance within the decade-long time period of this study. Only 2.3% of the areas 
burned by fires had subsequent MPB outbreaks. And the percentages were even smaller 
for fires and clearcuts in MPB infested forests, which were 1.7% and 0.8% of the area 
affected by MPB mortality, respectively. Although the probability of disturbance 
interactions was low in this region, their combined impact could potentially be different 
than the impacts of individual disturbances. As revealed from the ANOVA test, there was 
a significant difference in the magnitude of MPB and fire events between the fire+MPB 
target and population groups, but no significant difference was found between the 
MPB+fire target and population groups with regard to both events. Their comparison also 
revealed that stands experiencing fire+MPB interaction were attributed with higher fire 
magnitude but lower beetle outbreak magnitude than those without (Figure 9). 
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Figure 9. Differences in the magnitude of individual disturbances for areas with multiple 
disturbances. The boxes with solid color represent the target group, and the boxes with 
dash lines represent the random samples selected from the whole population besides the 
target group. The bottom and top of the box are the first and third quartiles, and the band 
inside the box is the median. The whiskers represent the 1.5 interquartile range of the 
lower and upper quartile. 

We outlined the major successional pathways in SRME (Figure 10). During the 13-year 
period of this study, the forest area recovering from all kinds of disturbances was 2,539 
km 2 , which accounted for 24.2% of the total disturbed area. To understand what was 
happening in the disturbed regions showing no sign of recovery till year 2011, we 
summed the area for each specific event according to their onset years (Figure 11), and 
found that most of those regions just experienced very recent disturbances. For instance, 
64.6% of the disturbances with no regrowth were happened within 5 years (Figure 11). 
This implies that the recovery rate may be higher than 24.2% since most regions were 
still undergoing disturbances. 
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Figure 10. Disturbance and recovery pathways in South Rocky Mountain Ecoregion. 
Percentages in the arrow indicate the percentage of area in the left box to that in the right 
box. 


For places that experienced a single disturbance, burned areas had the most regrowth, 
with a total area of 1,428 km 2 . Sites with only standing-replacing clearcut disturbances 
had the second largest area in regrowth (338 km 2 ), which was 49.8% of the total clearcut 
area, and MPB disturbed forests had the least amount of area classified as recovering 
(662 km 2 ). The regeneration rate, which was defined here as the ratio of post-disturbance 
regrowth area to the disturbed area, was 57.9% for fire burned sites, and 49.8% for 
clearcut sites over the 13-year study. Flowever, only 9% of the forest area was recovering 
within the MPB-induced mortality zones. In terms of multiple disturbances, their 
compounded effects on the subsequent recovery were evident. The results showed that 
the regeneration rate for fire+MPB was 14.9%, and MPB+fire produced the highest 
regeneration rate (64.9%) amongst all kinds of singular and compounded disturbances. 
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Figure 11. The summed disturbed area of mountain pine beetle (MPB), fire and clearcut 
in regions without regrowth according to their onset year. 

5 Discussion 

5.1 New classification framework and uncertainties 

Mapping over a decade of forest DR history at the ecoregion scale using moderate spatial 
resolution satellite imagery presented challenges that led to our development of new 
methodological techniques. Considering the large-scale geographic coverage, high levels 
of landscape heterogeneity and the diversity of disturbance agents, the mean OA of 85.9% 
demonstrated the efficiency of this framework. This framework can also contribute to the 
remote sensing community from two perspectives. One is the SOA/NOA sampling 
strategy. Selecting samples in the overlapping areas between adjacent Landsat scenes, in 
particular between paths on the same row, can reduce the time required for image 
interpretation or on-site data collection. By further including NOA clusters sharing 
dissimilar spectral features with the SOA clusters, we increased the overall spectral 
variability spanned by the selected samples and consequently enhanced the sample 
representativeness that is of particular importance for parametric algorithms (Gong et ah, 
2013). 

Another contribution exists in our usage of segments as the basic classification unit and 
the integration of spatio-temporal neighborhood features. Classification of segments can 
reduce predicted illogical annual changes since the labels are consistent over the course 
of the segment. The inclusion of an adequate number of spatio-temporal neighborhood 
features can improve the robustness of RF performance, but also minimize the “salt-and- 
pepper” effects from both spatial and temporal perspective. 
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Nevertheless, there are several limitations in our mapping procedure that might impact 
the subsequent ecological analysis. First, the temporal segmentation was run on a unified 
set of parameters that were calibrated with data from Grand County, Colorado (Liang et 
ah, 2014a), where a high concentration of both lodgepole pines and MPB mortality are 
located. The differences in MPB populations and the type and density of host tree species 
might impact the segmentation outputs, and result in a certain amount of false positive 
errors in MPB mortality pixels. Second, segments showing slow and gradual declining 
trends do not necessarily correspond to MPB-induced mortality. Forests experiencing 
drought-induced tree die-off, disturbance from other insects or disease-causing pathogens 
and mistletoes, or undergoing thinning and partial harvest, could also present this kind of 
pattern, and therefore be incorrectly attributed to MPB-mortality. Restricting our analysis 
within the MPB host tree species can minimize but not entirely avoid this confusion. 

5.2 Interactions between fire and MPB outbreaks 

Although the area of post-MPB fires was almost double that of post-MPB clearcuts, its 
coverage was fairly small compared with the vast area affected by MPB mortality. Thus, 
it remains challenging to draw any solid conclusions about the role of MPB mortality in 
increasing the risks of fire hazards in SRME. There is a growing body of ecological 
research exploring the relationships between MPB and fire, and some have conflicting 
conclusions (Hicke et ah, 2012; Schoennagel et ah, 2012; Simard et ah, 2011). Most of 
them were carried out by measuring the surface and canopy fuels in a time-since-beetle- 
outbreak chronosequence at a limited number of field sites to parameterize fire simulation 
models, and to predict future potential fire behavior. However, beetle-fire interactions 
could be nonlinear under different environments (Simard et ah, 2011), and the small 
quantities of on-site data used in past studies may not fully capture the underlying 
heterogeneity and variation. Our satellite measures, although lacking the details of the 
site-level collections, are capable of providing infonnation about where and when those 
events happened at the landscape level. 

In SRME, bum severity was largely unrelated to pre-fire MPB outbreak severity. For 
fires that burned through beetle-impacted forests, their severity did not show any 
significant difference from fires that burned through forests that were not affected by 
beetles. Beetles can affect the forest ecosystem in a variety of ways. Wind penetration 
through stands may increase if the beetle killed trees started to fall, but we did not 
observe a large proportion of fallen trees in this region during field visits and on high- 
resolution images. The surface and canopy fuel would be altered after the mortality of 
individual trees, but not in a consistent way. Shortly after beetle mortality, canopy fuels 
are high due to the declining foliage moisture. Over several years, the fuel load in the 
canopy layer may be reduced because limbs and trees fall, whereas the surface fuel loads 
would increase because of the thicken fuel bed and accumulated fallen dead needles. 

Thus, it is hard to predict the compound effects on fire potential since their contrary 
characteristics may be counteracting each other and leading to the uncorrelated 
relationship that we observed. Besides, fuel is not the only factor and sometimes there are 
other less important factors that control the burn severity. Studies using either satellite or 
field measures of bum severity indicated that relative humidity, wind speed, topography 
and burning conditions are often more critical factors in determining burn severity than 
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pre-fire outbreak severity (Bigler et al., 2005; Collins et al., 2007; Harvey et al, 2013; 
Romme et al., 2006). 

However, post-fire beetle severity has been shown to decrease with burn severity, which 
suggested that MPB usually cause less damage in forests in which there had been a prior 
fire event. Wildfire can impose immediate stress on trees. A prominent effect of fire 
injury is reflected in the low substrate quality, as the brood production declined with the 
degree of fire injury (Powell et al., 2012). Hence, heavily burned stands would harbor 
less density of beetles, and result in milder outbreaks. Moreover, fires consume a larger 
amount of biomass, and if the fire-affected stands did not recover back to their healthy 
status before the beetle attack, the decline in canopy cover will be less evident than when 
beetles attack more healthy stands. 

5.3 Impacts of singular and compounded disturbance on forest resilience mechanisms 

Our analysis found that more than half of the burned area was classified as regrowth. 
Lodgepole pines, which dominate our study area in the SRME, have serotinous cones, 
which need high temperatures to trigger their cones to open and release their seeds. This 
fire-dependent characteristic allows for post-fire tree regeneration soon after fires (Turner 
et al., 1997), although the abundance might vary with fire regimes and pre-fire seed 
density (Schoennagel et al., 2006). Lodegpole pine start seed production at an early age, 
approximately 5-10 years (Burns & Honkala, 1990), and thus it is not surprising to expect 
this high regeneration rate. Another question that many researchers and stakeholders are 
concerned with is how forest responds to post-MPB disturbances. Our results showed that 
a small portion of the MPB disturbed forests was on a trajectory to full recovery. The 
remaining areas were either under disturbance or were recovering at a slower rate. These 
observations may be related to dead trees dominating the spectral signal observed by 
Landsat, or to warmer and drier climate conditions that have been unfavorable for tree 
regeneration in recent years, as well as related to the size of disturbed patches being much 
larger than the short dispersal distances of tree species in the SRME (Bentz, 2010; Kurz 
et al., 2008; Sambaraju, 2012). 

Interactions between disturbance processes can amplify or mute the overall effects of 
changes (Vose et al., 2012), and the outcome is largely determined by the magnitude and 
direction of effects they created, as well as the order of events. One interesting piece of 
evidence we found is that beetle outbreaks could alter ecosystem response to subsequent 
fire by increasing the post-fire regeneration rate of serotinous lodgepole pines in SRME. 
Despite the fact that beetles might reduce seed availability by killing large seed- 
producing trees (Bjorklund & lindgren, 2009), they can also help in redistributing the 
remaining viable cones and seeds, like what the wind blow or animal carrier do 
(Chambers & MacMahon, 1994), and promote the generation of a larger number of 
seedlings as aided by the following fire. However, this finding could only be valid in this 
particular ecosystem and under specific conditions, such as stand age, density and bum 
severity, because we found some contradictory conclusions in other studies. For instance, 
Harvey et al. (2014) suggested that there is no direct effect of outbreak severity on initial 
post-fire regeneration of lodgepole pine in Greater Yellowstone. 

6 Conclusions 
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In this study, we documented the disturbance and recovery history in the lodgepole pine 
forests in the SRME using Landsat time series classification maps, which were generated 
from a novel SOA/NOA sampling strategy, the integration of Landtrendr with an RF 
classifier to attribute disturbance type, and the inclusion of spatio-temporal neighbor 
features. By tracking the spatial and temporal of patterns of MPB mortality and burned 
areas, we found that the linkage, strength, and direction of interactions between bark 
beetle mortality and fire to be related to their sequence of occurrence. The severity of fire 
was unlinked to pre-fire outbreak severity, and no evidence was found towards the role of 
MPB in increasing the risks of fire hazard. However, post-fire beetle severity tended to 
decrease in more severe fire injured stands. Most recovery events happened immediately 
after stand-replacing disturbances, whereas chronic beetle disturbance that contributed 
most to the forest loss, showed an extremely low regeneration establishment rate. 
Although the occurrence probability of multiple disturbance events is low in this region, 
their compound effects on the resilience mechanism of serotinous lodgepole pine forests 
can be of great ecological value. The most evident impact is the increased regeneration 
rate after the combination of pre-fire beetle outbreaks and subsequent fire. As bark beetle 
outbreaks and wildfire occurrence are both predicted to increase with continued climate 
wanning in North America (Kurz et ah, 2008; Westerling et al. 2011), the maps and 
findings presented in this study could serve as baseline information for assessing 
landscape-scale effects of disturbance and recovery on forest ecosystems, and help 
infonn effective management options. 
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Chapter 4 Projecting future potential patterns of beetle outbreaks in the 
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Abstract 

The recent widespread mountain pine beetle (MPB) outbreak in the Southern Rocky 
Mountains presents an opportunity to investigate the relative influence of anthropogenic, 
biologic, and physical drivers that have shaped the spatiotemporal patterns of the 
outbreak. The aim of this study was to quantify the landscape-level drivers that explained 
the dynamic patterns of MPB mortality, and simulate areas with future potential MPB 
mortality under projected climate-change scenarios in Grand County, Colorado, USA. 
The outbreak patterns of MPB were characterized by analysis of a decade-long Landsat 
time-series stack, aided by automatic attribution of change detected by the Landsat-based 
Detection of Trends in Disturbance and Recovery algorithm (LandTrendr). The annual 
area of new MPB mortality was then related to a suite of anthropogenic, biologic, and 
physical predictor variables under a general linear model (GLM) framework. Data from 
years 2001-2005 were used to train the model and data from years 2006-2011 were 
retained for validation. After stepwise removal of non-significant predictors, the 
remaining predictors in the GLM indicated that neighborhood mortality, winter mean 
temperature anomaly, and residential housing density were positively associated with 
MPB mortality, whereas summer precipitation was negatively related. The final model 
had an average area under the curve (AUC) of a receiver operating characteristic plot 
value of 0.72 in predicting the annual area of new mortality for the independent 
validation years, and the mean deviation from the base maps in the MPB mortality areal 
estimates was around 5%. The extent of MPB mortality will likely expand under two 
climate-change scenarios (RCP 4.5 and 8.5) in Grand County, which implies that the 
impacts of MPB outbreaks on vegetation composition and structure, and ecosystem 
functioning are likely to increase in the future. 

1 Introduction 

As a native species, mountain pine beetle ( Dendroctonus ponderosae ; MPB) populations 
have existed at endemic levels and periodically have grown to epidemic levels in the pine 
forests of western North America for centuries (Amman, 1977, Baker & Veblen, 1990; 
Raffa et ah, 2008). By infesting and killing older and stressed trees with larger diameters, 
MPB plays a critical role in shaping forest composition and structure, accelerating the 
movement of nutrients in biogeochemical cycles, and affecting forest productivity 
(Collins etal., 2011;Edburg et ah, 2012). In recent decades this historical balance has 
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been disrupted, and the area affected by MPB has been vastly extended, exceeding the 
extent and impacts of outbreaks documented in the past 125 years (Raffa et ah, 2008). 
The current MPB outbreak has impacted large expanses of lodgepole and ponderosa pine 
forests, reduced their ability to act as carbon sinks (Caldwell etal., 2013, Kurz et ah, 
2008; Running, 2008), altered wildfire hazards (Hicke etal., 2012a, Jenkins etal., 
2008, Parker etal., 2006; Schoennagel etal., 2012), modified local surface energy 
balance (Boon, 2009), threatened water quality (Mikkelson etal., 2013), and changed 
regional climate (Maness et ah, 2013). 

The population dynamics of bark beetles are governed by a variety of biotic and abiotic 
factors and their interactions (Raffa et ah, 2008). Forest characteristics, including 
homogenous even-aged, high-density, and large-diameter stands, are favorable for MPB 
mass attack (Raffa & Berryman, 1983). Long-term drought or other events causing stress 
can exert either positive or negative effects on tree susceptibility to beetle attack, and the 
overall impact remains controversial. While the primary defense mechanism of trees will 
be weakened by drought stress because of reduced resin quantities (Creeden et al., 
2014, Preisler etal., 2012; Raffa etal., 2008), beetle brood production can also be 
reduced since the tree's phloem thickness is attenuated (Amman & Cole, 1983). Thennal 
regimes, typically represented by minimum winter temperature or year-round 
temperature, impact beetles' developmental timing, cold-induced mortality, and the 
associated fungal community (Bentz etal., 2010; Preisler etal., 2012). Meanwhile, 
factors like elevation, direct solar radiation, and beetle mortality in adjacent areas have 
also been indicated as important predictors of outbreaks (Coops et al., 2006b, Simard 
etal., 2012, Walter & Platt, 2013;Wulder etal., 2006b). Most previous studies have 
focused on the effects of one or several factors. In this study, we took a comprehensive 
approach and considered a large set of relevant factors to improve our understanding of 
the spatiotemporal patterns of MPB outbreaks and investigate their drivers. 

For this study, we were also concerned about the identification of forested areas with 
high risk for future MPB mortality. There have been a number of prior efforts to predict 
patterns of MPB mortality. Through an integrated seasonality and cold tolerance 
model, Bentz etal. (2010) suggested that rising temperatures could increase MPB 
population growth rates, and their range would expand along latitude and elevation 
gradients. Aided by an ecological niche model, Evangelista et al. (2011) predicted that 
new areas of forest susceptible to MPB mortality would emerge over time but the existing 
area of susceptible forests to MPB mortality would also shrink, leading to an overall 
decrease in the amount of suitable habitat area for MPB in the future. Using a process- 
based model, Hicke et al. (2006) found that projected warming in the western United 
States will result in substantial reductions in the overall area of adaptive seasonality (the 
synchronous emergence of adults that allows MPB to overwhelm tree defenses). Unlike 
population models, which can improve the mechanistic understandings of biological 
responses to environmental variability, but may consider a limited number of explanatory 
variables because of model complexity, statistical approaches are capable of 
incorporating a large number of explanatory variables and quantifying their relative roles. 
This is a crucial preliminary step before adopting and improving process-based 
mechanistic models. 
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Understanding the factors driving patterns of MPB outbreaks and predicting future 
outbreaks has been challenging given the types of data available that depict the spatial 
and temporal extents of outbreaks. The quality of response variables could affect the 
performance of predictive models. In general, locations of MPB mortality are collected in 
the field or extracted from remotely sensed imagery. Although in-situ surveys can 
provide accurate data, they often have restricted geographic and temporal coverage. State 
and federal agencies have conducted Forest Health Monitoring Aerial Detection Surveys 
(ADS) to identify forest disturbances since the mid-twentieth century (Man, 2010). These 
publicly available datasets have been used extensively in many fields, but errors 
introduced via observer fatigue, observer-to-observer variation, misregistration and the 
scale of observation (Meigs et ah, 2011) are rarely estimated and could introduce an 
unknown amount of uncertainty. 

Remote sensing of forest disturbances offers an alternative to the ADS data for 
monitoring tree mortality caused by insect outbreaks (Coops et ah, 2006a; White et ah, 
2004). Landsat data are especially popular for this application because they are freely 
available, and have multispectral data, a broad spatial extent, and temporal continuity. 
For these reasons, Landsat time series stacks (LTSS) have been used in large-scale efforts 
to detect forest disturbances (Masek et ah, 2013) using the Vegetation Change Tracker 
algorithm (VCT; Huang etal., 2010). Remotely sensed disturbance maps produced by 
VCT and similar change-detection algorithms like the Landsat-based Detection of Trends 
in Disturbance and Recovery algorithm (LandTrendr; Kennedy et ah, 2010) currently 
lack infonnation about the cause of the disturbance, and that limits their utility for use in 
our study. Notwithstanding, several studies have demonstrated the usefulness of Landsat 
in capturing the patterns of MPB-caused tree mortality at various geographic scales 
(e.g. Masek et ah, 2013; Meddens et ah, 2013). Considering the uncertainties in the ADS 
data, and limitations of existing VCT and LandTrendr change-detection products, we 
utilized data from an automated procedure that labeled disturbance types (especially 
MPB mortality) detected by LandTrendr in an LTSS to generate spatially explicit annual 
maps of MPB occurrences over a decade-long time span (Liang et ah, 2014a). 

In this paper, we integrated remote sensing techniques and statistical models to evaluate 
the effects of a set of factors affecting the dynamic pattern of MPB mortality, and 
projected future MPB mortality in response to climate change. Our aim was to address 
the following questions: What drivers promote the extensive development and 
progressive MPB outbreak in an area situated in the Southern Rocky Mountains 
ecoregion? How accurately can we predict MPB disturbance with this set of response and 
explanatory variables? And what will future outbreak trends be? 

2 Study area 

Grand County is located in north central Colorado, covering approximately 4830 square 
kilometers of the Southern Rocky Mountains ecoregion (Figure 1). The elevation ranges 
from 2225 m along the Colorado River to 4131 m at the summit of the Continental 
Divide (Grand County Department of Natural Resources, 2006). Its climate is 
characterized by year-round sunny days (around 244 days/year on average), with average 
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summer temperatures of 26.6 °C, and the average rainfall of approximately 30.48 cm 
(Grand County Department of Natural Resources, 2006). The diversity of elevation, soil, 
climate, as well as strong topographic-moisture gradients leads to a variety of vegetation 
composition within the county, among which sagebrush shrub and steppe is the most 
dominant ecosystem. Lodgepole pine forests occupy a quarter of the landmass, followed 
by spruce-fir forests and aspen forests (Grand County Department of Natural Resources, 
2006). In recent decades, MPB infestation, wildfire, and timber harvesting are recognized 
as the three major disturbance agents in Grand County. Wildfire occurrence has been 
low, but the widespread MPB outbreak affected approximately 68% of privately owned 
land and 70% of federally owned land (Witcosky, 2007). 
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Figure 1. Spatial location of (a) the State of Colorado within the United States of 
America, (b) Grand County within Colorado, and (c) areas of lodgepole pine forest 
(derived from the LANDFIRE existing vegetation type data layer) within Grand County, 
Colorado. 

3 Methods 

3.1 Change detection analysis in detecting long-term MPB outbreaks 

Maps of MPB mortality in Grand County were generated by automatic attribution of 
LandTrendr segmentation outputs applied to a time series of 17 Landsat images spanning 
2000-2011 (path 34, row 32; Liang et ah, 2014a). This approach integrated a temporal 
segmentation technique to identify areas with change (Kennedy et ah, 2010) and a 
decision tree modeling procedure to attribute the changes to specific disturbance types 
(Liang et ah, 2014a). The steps in this approach were to enhance the signal-to-noise ratio 
of the time series curve via linear regressions, and decompose it into a sequence of 
straight-line segments, whose event attributes were identified based on the segment 
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characteristics (duration, magnitude and vertex value) using calibrated decision tree rules 
(Figure 2). The temporal trajectories were constructed using the Nonnalized Bum Ratio 
(Key & Benson, 2006). Key parameters that were used to calibrate both temporal 
segmentation and decision tree components were trained using ground-truth data from 
106 sample locations, visually selected and interpreted from 1-m resolution U.S. 
Department of Agriculture (TJSDA) National Agricultural Imagery Program (NAIP) 
imagery (available in 2005, 2009, and 2011 from USDA Geospatial Data Gateway). 
Technical details can be found in Liang et al. (2014a). The trend analysis outputs 
included a series of disturbance maps that depict healthy forest, forest with MPB 
mortality, and clearcut areas at an annual time step. The disturbance maps were validated 
with a set of randomly placed NAIP test samples in Liang et al. (2014a). 



Figure 2. Processing steps in the change-detection analysis. LTSS: Landsat time-series 
stacks; NAIP: National Agricultural Imagery Program; LEDAPS: Landsat Ecosystem 
Disturbance Adaptive Processing System; FMASK: Function of Mask; LandTrendr: 
Landsat-based Detection of Trends in Disturbance and Recovery algorithm; NBR: 
Normalized Burn Ratio; MPB: mountain pine beetle 

3.2 Model development 

3.2.1 Response variable 

We were interested in simulating the spread of MPB mortality into new regions over 
time, instead of modeling their suitable habitats. Thus, we set the response variable in our 
models to be the annual presence (case) of new MPB mortality and absence (control) of 
MPB mortality, where presence refers to the new mortality in a pixel which was healthy 
in the previous year. The time-series of disturbance maps were used as the base for 
sample selection. First, we made a random sample of newly emerged areas of MPB 
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mortality stratified by each year from 2001 to 2011. A stratum containing persistently 
(2001-2011) healthy forest pixels was also constructed for control sample selection. 
Second, since the sample size in each stratum detennines the distance between 
observations, and thus affects the potential for spatial autocorrelation to influence model 
interpretation, we initially selected 300 sample units from each stratum, and successively 
decreased the sample size with decrements of 10. For each sample size, a general linear 
model (GLM) was fit and Moran’s / was calculated to test for spatial autocorrelation in 
the model residuals (Moran, 1950). We ultimately selected the largest sample size that 
had insignificant spatial autocorrelation in the model residuals, that is, we selected a 
sample with enough spacing between points so that spatial autocorrelation effects were 
avoided. 

A set of training samples was extracted using observations from 2001 to 2005, and the 
model with forward predictions was validated with observations from 2006 to 2011. Year 
2000 was not included in modeling because we could not infer its prior infonnation for 
certain predictors, such as distance to nearest cell with MPB mortality in the previous 
year. Validation data collected after 2008, when spread of the MPB outbreak was 
reduced, were used to determine whether or not the model over predicted MPB mortality. 

3.2.2 Explanatory variables 

Thirty-four biotic and abiotic variables that fell into seven categories were used to 
develop a number of GLMs (Table 1). All datasets, except for the two climate variables, 
were in raster fonnat with the same spatial resolution as our disturbance maps (30 m). 

The values from all the rasters were extracted at the locations of the training or validation 
sample points. 

Anthropogenic variables included residential housing density and distance to nearest road. 
Both of them are proxies for the intensity of human activities, whose impacts have rarely 
been investigated in prior MPB associated studies. The biggest human impact on forest 
ecosystems is likely to be habitat fragmentation since silvicultural treatments, such as 
thinning or logging, are still the most common management strategy in mitigating MPB 
outbreaks (Coops et ah, 2008). How human intervention affects MPB host selection 
remains unknown, however. To quantify the magnitude of this potential effect, we used 
residential housing density data from the 2010 U.S. Census Bureau block-level housing- 
density data (Radeloff et ah, 2010), and distance to nearest road from the National 
Overview Road Metric Euclidean Distance dataset (Watts et ah, 2007). Both continuous 
variables were log + 1 transfonned prior to use because they had skewed distributions. 

Lodgepole pine forests exist along a topographic-moisture gradient that controls 
vegetation growth as a function of soil water holding capacity, evapotranspiration and 
surface runoff. Six topographic variables derived from the Shuttle Radar Topography 
Mission Digital Elevation Model were used to represent this gradient (USGS, 2004). 
Aspect was recalculated in a way that the first 45° from true north to east was recorded as 
1, and increased by 1 as the aspect increased every 45° clockwise. Southwestness is a 
cosine-transformation of aspect that ranges from -1 to 1 (Franklin et ah, 2000). 
Topographic wetness index (TWI) is the steady-state humidity index (Beven & Kirkby, 
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1979). Higher TWI values indicate greater soil moisture, and we expected TWI and the 
distance to stream channel networks to be negatively correlated with MPB mortality. 

We used tree cover prior to the MPB outbreak as a proxy for the pre-disturbance health 
and abundance of MPB host species. The Landsat vegetation continuous fields tree cover 
layer for circa-2000 provides estimates of the aboveground woody vegetation percentage 
in each 30 m pixel (Sexton etal., 2013). Under most circumstances, canopy cover is 
positively and significantly correlated with diameter at breast height (Gill et ah, 2000), 
and thus high canopy cover usually represents large-diameter trees, which are more likely 
to be attacked by MPB (Amman, 1977; Klutsch et ah, 2009). 

Spatial proximity to areas previously affected by MPB can be a critical facilitator in 
driving local outbreaks because MPB is a relatively poor disperser (Simard et ah, 2012), 
and we included two types of dispersal-related variables to represent it. Distance to 
nearest mortality is computed as the Euclidean distance from one cell to its closest cell 
with MPB mortality in the previous year. Shorter distance between sites enhances their 
connectivity, and thus increases the probability of beetle dispersal to adjacent healthy 
sites. 

Neighborhood mortality is the amount of adjacent pixels with MPB mortality in the 
previous year. Because spatial synchrony is prevalent among beetle populations during 
epidemic years (Aukema et ah, 2006), more beetle presence in the immediate 
neighborhood increases the likelihood of a mass attack on adjacent healthy forest. The 
scale of neighborhood depends on the beetles' dispersal pattern, which has been 
summarized into two modes: short-distance and long-distance dispersal. Short-distance 
dispersal happens within stands (Safranyik et ah, 1989), and long-distance dispersal 
usually occurs when beetles are transported above the forest canopy by wind (de la 
Giroday et ah, 2011; Jackson et ah, 2008; Robertson et ah, 2007; Robertson et ah, 2009). 
The common distances in the short-distance range dispersal are 30-50 m (Robertson 
etal., 2007, Safranyik etal., 1989; Safranyik etal., 1992), whereas the long-distance 
flight dispersal that depends on the wind speed, preflight weight, flight duration, and lipid 
content (Evenden et ah, 2014) can be more variable, ranging from several to tens of 
kilometers. In field observations, 2-3 km were commonly found to be the maximum 
distance beetles can disperse by entering a new stand from surrounding areas (Robertson 
et ah, 2008; Robertson et ah, 2009), whereas laboratory flight mill bioassay showed that 
the mean MPB flight distance ranged between 2.12 and 5.95 km (Evenden et ah, 2014). 
Since there is no consensus about which mode is more important in driving the beetle 
expansion, we defined a number of neighborhood distances: 30 m, 100 m, 250 m, 500 m, 
1 km, 1.5 km, 2 km and 3 km. Those distances were used as the radius of a circular 
window, and all pixels with MPB mortality in the previous year covered by this window 
would then be counted to be the neighborhood mortality value for the center cell. 
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Table 1. Predictor variables selected for use in the general linear models (GLMs). 


Category 

Variables 

Abbreviation 

Resolution 

Unit 

1. Anthropogenic 

Distance to the nearest road 

road 

30 m 

m 


Residential housing density 

house 

30 m 

m 

2. Topography 

Elevation 

dem 

30 m 

m 


Aspect 

aspect 


degree 


Slope 

slope 


degree 


Southwestness 

sw 


no unit 


Topographic Wetness Index 

twi 


no unit 


Distance to channel network 

dis2chan 


m 

3. Vegetation condition 

Tree cover 

tc 

30 m 

percentage 

4. Spatial proximity 

Distance to closest mortality in the previous year 

dis2prev 

30 m 

m 

5.Neighborhood 
mortality 

Number of pixels in an 8-pixel neighborhood 
with MPB mortality in the previous year 

Nm_30m 


no unit 


Number of pixels in a 100 m circular 
neighborhood with MPB mortality in the 
previous year 

nm 100m 


no unit 


Same as above, but in a 250 m circular 
neighborhood 

nm 250m 


no unit 


Same as above, but in a 500 m circular 
neighborhood 

nm 500m 


no unit 



Same as above, but in a 1 km circular 
neighborhood 

Same as above, but in a 1500 m circular 
neighborhood 

Same as above, but in a 2 km circular 
neighborhood 

Same as above, but in a 3 km circular 
neighborhood 

6. Climate Summer mean temperature* 

Winter mean temperature* 

Summer mean precipitation 

Wannest temperature 

Coldest temperature from October to May 

Mean annual temperature of previous year 

Mean annual precipitation of previous year 
Mean summer precipitation of previous year 

7. Climate anomalies Mean summer precipitation anomaly 

Mean winter temperature anomaly 


O') 

00 


nm 1km 


no unit 


nm 1500m 


no unit 

nm 2km 


no unit 

nm 3 km 


no unit 

tmean summer 

4 km 

°C 

tmean winter 


°C 

ppt_summer_cur 


°c 

Tmax cur 


°c 

Tmin cur 


°c 

tmean last 


°c 

ppt mean last 


mm 

ppt summer last 


mm 

ppt_summer_cur2nor 

mal 

4 km 

mm 

tmean winter2normal 


°C 



Wannest temperature anomaly 

tmax cur2nonnal 

°C 

Coldest temperature anomaly 

tmin cur2normal 

°C 

Mean annual precipitation anomaly of previous 
year 

ppt mean last2norma 

1 

mm 

Mean annual temperature anomaly of previous 
year 

tmean last2normal 

°C 

Mean summer temperature anomaly of previous 
year 

tmean summer2nonn 
al 

°C 

Mean summer precipitation anomaly of previous 
year 

ppt summer last2nor 
mal 

mm 


* Summer is from June through August, and winter is from December through February. 



We used climate datasets generated by the Parameter-evaluation Regressions on 
Independent Slopes Model (PRISM Climate Group, 2010) in which monthly and annual 
weather data are available at a resolution of approximately 4 km (Daly et ah, 2002). 
Besides annual mean temperature and precipitation infonnation in the PRISM data, we 
derived six additional variables that are controlling factors in the beetles' life cycle 
(Kaufmann et al., 2008): summer mean temperature, winter mean temperature, summer 
mean precipitation, warmest summer temperature, coldest winter temperature, and 
summer mean precipitation in the previous year. Additionally, we computed eight climate 
anomalies by taking the differences between each climate variable and their 30-year 
averages from 1981 to 2010, since climate change has been indicated to have both direct 
and indirect impact on MPB outbreaks (Kurz et al., 2008). 

3.2.3 Modeling approach 

We used general linear models (GLM) with a logit link and binary response to identify 
which variables explained recent patterns of MPB mortality and assess potential new 
areas of MPB mortality. We first applied univariate GLMs to each of the 28 predictor 
variables to assess their individual relationship with MPB mortality. These models were 
evaluated by their coefficient estimates and associated significance tests. Their spatial 
autocorrelation effects were examined using Moran's/on the model residuals. Variables 
with /7-values greater than 0.05 on the spatial autocorrelation tests were excluded from 
further analyses, to avoid violation of the assumption of independently and identically 
distributed errors in GLMs (Donnann et al., 2007). We then fit a full model starting with 
all predictor variables that were not removed because of spatial autocorrelation in the 
univariate models. We retained statistically significant predictor variables in the full 
model based on the Bayesian Information Criterion (BIC) in a multiple backward 
stepwise selection algorithm implemented in R (R Core Team, 2013). We chose BIC 
because of the large number of predictor variables and BIC penalizes the number of 
parameters more strongly than commonly used Akaike Information Criterion. 


Predictive maps of MPB mortality were generated using the final GLM after backward 
stepwise selection with the equation: 


P ~ 1/(1 + eX P(~(Po + Pl%l + P 2%2 3 -b PiX t + C))) 


where P is the probability of MPB mortality, (B 0 is the intercept, Xps a predictor variable, 
and Pi is the regression coefficient estimate for the associated predictor variable X;. The 
correction factor (C) was used to account for model bias introduced because of different 
ratios of case and control observations in the sample and in the population (Manly et al., 
2002 ): 


C = log ( 


# sample, control / # pop. controh 
# sample, case / # pop. case > 


where # sample.control and # pop.control are the number of healthy (no MPB mortality) 
occurrences in the sample and the population, respectively; # sample.case and # pop.case 
are the number of observations with MPB mortality in the sample and the population, 
respectively. 
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3.3 Model performance evaluation 

Model performance was evaluated in three ways. (1) The new presence of MPB mortality 
each year predicted by our GLM was visually compared with observed MPB mortality in 
the Landsat classifications. This provided an immediate and intuitive way to perfonn the 
evaluation. (2) The areal estimates of MPB mortality derived from remote sensing and 
GLM were compared and their accuracy was evaluated against the criteria: the peak of 
tree mortality should occur around 2005 and 2008 (Klutsch et ah, 2009), and the rate of 
change should be the greatest at the beginning and then reduced until reaching a stable 
level. (3) A quantitative accuracy assessment was also conducted. First, 10-fold cross- 
validation was applied to test the model performance in the training phases to avoid 
problems like overfitting. Second, to detennine the model's predictive capacity, predicted 
MPB mortality was compared with an independent dataset extracted from the base maps 
from 2006 to 2011 in two ways. A set of 3000 points were randomly picked for areas of 
healthy forest and MPB mortality from the base image of each year, in order to provide 
an overall evaluation of landscape-level patterns. Another set of 3000 points were 
selected from areas of persistent healthy forest and new MPB mortality each year, to 
access the prediction ability for the occurrence of new MPB mortality. Both evaluations 
were judged by the area under curve (AUC) of the receiver operating characteristics 
(ROC) curve, as well as overall accuracy (OA; Fielding & Bell, 1997; Hanley & McNeil, 
1982). Overall accuracy was generated from the confusion matrix of the binary maps, 
where the probability of MPB mortality was separated from healthy forest with an 
optimal threshold calculated according to a criterion that maximizes the sum of 
sensitivity and specificity (Freeman & Moisen, 2008). 

3.4 Projections under future climate scenarios 

Future MPB mortality was projected using future climate conditions from the Coupled 
Model Intercomparison Project 5 (CMIP5). We used results from 14 global climate 
models (GCMs) downscaled to 4-km resolution by the Multivariate Adaptive Constructed 
Analogs (MACA) method, which was designed for wildfire applications in the western 
USA (Abatzoglou & Brown, 2012; Table 2). Projections under two future Representative 
Concentration Pathways, RCP4.5 and RCP8.5, have been adapted by MACA for 
downscaling. They represent a high pathway for which relative radiative forcing reaches 
>8.5 W/m 2 by 2100, and an intennediate pathway where radiative forcing is stabilized at 

4.5 W/m 2 after 2100, separately. We made projections for each GCM and RCP, generated 
spatially-explicit maps, and calculated mean, 25th quartile, 75th quartile, lowest and 
highest probability of all GCM projections for each year and each RCP. Because data 
incorporating the climate change effects on the future forest extent and residential 
housing density at the spatial resolution of our model are not currently available, our 
projections assumed no changes in the distributions of forest and housing densities, as 
well as other controlling factors, such as the supply of nutritionally optimal host trees 
(Raffa et al., 2008). 


Table 2. The 14 global climate models (GCMs) from which downscaled climate 
projections were used in this paper and the modeling groups developing them. 


Model Name 


Modeling Group 
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BCC-CSM1.1 

Beijing Climate Center, China Meteorological 
Administration 

BNU-ESM 

College of Global Change and Earth System Science, 
Beijing Normal University 

CanESM2 

Canadian Centre for Climate Modelling and Analysis 

CNRM-CM5 

Centre National de Recherches Meteorologiques / Centre 
Europeen de Recherche et Formation Avancee en Calcul 
Scientifique 

CSIRO-Mk3.6.0 

Commonwealth Scientific and Industrial Research 
Organization in collaboration with Queensland Climate 
Change Centre of Excellence 

GFDL-ESM2G 

GFDL-ESM2M 

National Oceanic and Atmospheric Administration 
Geophysical Fluid Dynamics Laboratory 

HadGEM2-ES 

HadGEM2-CC 

Met Office Hadley Centre (additional HadGEM2-ES 
realizations contributed by Instituto Nacional de Pesquisas 

Espaciais) 

INM-CM4 

Institute for Numerical Mathematics 

MIROC5 

Atmosphere and Ocean Research Institute (The University of 
Tokyo), National Institute for Environmental Studies, and 
Japan Agency for Marine-Earth Science and Technology 

MIROC-ESM 

MIROC-ESM- 

CHEM 

Japan Agency for Marine-Earth Science and Technology, 
Atmosphere and Ocean Research Institute (The University of 
Tokyo), and National Institute for Environmental Studies 

MRI-CGCM3 

Meteorological Research Institute 


4 Results 

In Liang et al. (2014a), we reported that our Landsat-based change-detection analysis for 
mapping disturbances resulted in an overall accuracy (OA) ranging from 87% to 94%, 
which was 20-30% higher than single-scene classifications performed by a maximum 
likelihood classifier and an ensemble random forest classifier. Because of Landsat's 
medium resolution, the percentage of dead trees in one pixel can be variable. We visually 
interpreted the dead tree cover percentage on NAIP imagery within the 30 m pixel 
window of the test samples, and found that 90% of the MPB-mortality pixels had more 
than 50% dead tree cover, while 58% of them have more than 80% dead tree cover. Since 
the Landsat-derived MPB mortality data were used as observations when constructing 
our models, our results should be interpreted as explaining and predicting the 
spatiotemporal patterns of moderate to severe MPB mortality (>50% dead tree cover). 
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Sixteen out of 34 predictors had statistically significant coefficients at the 0.05 level as 
tested by univariate GLMs (Table 3), and none of the predictors produced models with 
spatially autocorrelated residuals except TWI. Residential housing density, all 
neighborhood mortality variables, mean temperature in the previous year, mean summer 
precipitation and maximum summer temperature, and winter mean temperature anomaly 
were positively related with MPB mortality. Negative relationships were found between 
MPB mortality and elevation, distance to MPB mortality in previous year, mean annual 
precipitation in the previous year, and mean summer precipitation. 


Table 3. Predictor variable coefficients and significance levels, and Moran’s I for 
univariate general linear models (GLMs), and predictor variable coefficients and 
significance levels, and standard errors for the full multivariate GLM and the final 
multivariate GLM after stepwise selection. See table 1 for explanations of the predictor 
variable abbreviations. 


Variable 

Univariate GLMs 

Multivariate GLM 
(Ml) 

Multivariate 
GLM (final) 

abbreviation 






coef 

Moran’s I 

coef 

SE 

coef SE 

road 

0.126 

-0.026 

0.396** 

0.143 


house 

1.089* 

-0.034 

1.258 

0.661 

1.198* 0.477 

dem 

-0.003** 

-0.042 

0.002 

0.002 


aspect 

0.041 

-0.017 

0.038 

0.078 


slope 

-0.026 

-0.024 

-0.026 

0.025 


Sw 

0.136 

-0.017 

0.140 

0.264 


Twi 

1.104*** 

0.132* 




dis2chan 

0.001 

-0.019 

0.000 

0.004 


tc 

-0.020 

-0.025 

-0.009 

0.016 


dis2prev 

0.211*** 

-0.062 

-0.070** 

0.063 


nrn 30m 

0.582*** 

-0.028 

0.660 

0.205 

0.642*** 0.097 

nrn 100m 

0.085*** 

-0.045 

0.026 

0.047 


nrn 250m 

0.010*** 

-0.053 

-0.006 

0.010 


nrn 500m 

0.002*** 

-0.050 

0.001 

0.003 


nrn 1km 

0.001** 

-0.046 

0.000 

0.001 


nrn 1500m 

0.000* 

-0.041 

0.000 

0.001 
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nrn 2km 

0.000* 

-0.039 

0.002 

0.002 



nrn 3 km 

0.000* 

-0.039 

-0.003 

0.003 



ppt mean last 

-0.003** 

-0.040 

0.000 

0.004 



tmean last 

0.405** 

-0.046 

11.880 

7.005 



tmean winter 

0.132 

-0.013 

-5.131 

3.280 



tmean summer 

0.246* 

-0.031 

-6.126* 

3.619 



pptsummercur 

-0.014 

-0.026 

3.338* 

36.820 

-0.037** 

0.014 

pptsummerlast 

-0.011 

-0.029 

-3.395 

36.820 



tmin cur 

-0.014 

-0.021 

0.048 

0.786 



tmax cur 

0.155** 

-0.026 

-0.384 

0.870 



ppt mean last2normal 

0.000 

-0.020 

-0.018 

0.010 



tmean last2normal 

-0.529 

-0.012 

-13.10 

7.173 



tmean winter2normal 

0.030 

-0.018 

6.801 

3.404 

0.659** 

0.205 

tmean summer2nonnal 

-0.043 

-0.021 

8.662 

3.800 



ppt summer cur2normal 

0.002 

-0.020 

-3.388 

36.820 



ppt summer last2normal 

0.005 

-0.019 

3.459 

36.820 



tmin cur2normal 

0.033 

-0.017 

-0.313 

0.794 



tmax cur2nonnal 

-0.005 

-0.020 

-1.575 

1.087 



(intercept) 



15.330 

19.340 

0.460 

0.652 


Note: coef - coefficient estimate on the variable; SE - standard error; * denotes 
significance level of 0.05; ** denotes significance level of 0.01; and *** denotes 
significance level of 0.001. 

Four predictor variables were retained in the full model after the backward stepwise 
selection (Table 3). Residential housing density, number of pixels in the nearest eight 
pixels that had MPB mortality in the previous year, and winter mean temperature 
anomaly were positively related to the likelihood of MPB mortality. Summer 
precipitation was also retained in the full model, but was negatively related to MPB 
mortality. 

Besides quantifying the relative influence of the various drivers of MPB mortality, we 
were also curious about the predictive capacity of our GLM in a spatially explicit context. 
By assuming the satellite-derived disturbance maps were a true representation of 
landscape patterns of MPB mortality, we compared them against the GLM predictions for 
the validation years (Figure 3). We observed that the predicted areas of MPB mortality 
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generally matched well with the Landsat-based observations across the landscape. For 
instance, the northeast corner of Grand County is an area where the MPB outbreak 
progressively grew from 2006 to 2011. By carefully examining this zone, we found that 
the shrinking extent of healthy forests and the spread of MPB mortality as predicted by 
our model were basically in accordance with the satellite observed patterns. In the 
meantime, the areal estimates of MPB mortality predicted by the GLM for the whole 
county deviated little from the Landsat base maps, with 13% as the largest relative 
difference and the smallest relative difference was only 2%. Only year 2006 was under 
predicted and the remaining five independent validation years were over predicted. The 
annual predictions followed a pattern similar to that of the observed data (Figure 4): both 
showed a steady increase in the amount of area with MPB mortality, with a more rapid 
rate of increase before 2008 and slower rate of change after that. In terms of quantitative 
evaluations, the average AUC generated from 10-fold cross-validation was 0.97, with 
small variation among years. Overall accuracy ranged from 82 to 93% and the average 
was 88% (Table 4). When we assessed accuracy only in areas of new MPB mortality 
each year, our model achieved a mean AUC of 0.72, a mean OA of 0.66. 



2006 2007 2008 2009 2010 2011 



Figure 3. The comparison between observed disturbances in the change-detection 
analysis (panels a, c) and predictions of MPB mortality using the final GLM (panels b, d) 
in the independent validation years. Panels c and d show the detailed images from the 
northeast corner of the study area, the area of which is indicated by the circle in the 
leftmost panel a. 
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Figure 4. Accumulated area of mountain pine beetle (MPB) mortality detected by 
Landsat and predicted by our general linear model from 2001 to 2011, as well as 
forecasted MPB mortality to 2050. 


Table 4. Comparison of the final general linear model (GLM) to Landsat observed 
disturbances. _ 

AA1 AA2 




AUC 

OA 

AUC 

OA 

Landsat 

GLM 

RD 







(km 2 ) 

(km 2 ) 


Training 

period 

2000- 

2005 

0.98 

0.87 



833.23 

861.62 

0.03 


2006 

0.94 

0.89 

0.67 

0.63 

950.03 

934.55 

-0.02 

Validation 

2007 

0.96 

0.87 

0.63 

0.58 

1031.18 

1079.10 

0.05 

years 

2008 

0.97 

0.91 

0.66 

0.57 

1077.13 

1096.70 

0.02 


2009 

0.98 

0.82 

0.77 

0.73 

1087.08 

1230.60 

0.13 


2010 

0.98 

0.93 

0.79 

0.73 

1088.95 

1126.70 

0.03 


2011 

0.98 

0.87 

0.80 

0.71 

1082.05 

1168.40 

0.08 


AA1: accuracy assessment on the predicted annual image (previous year's + new 
mortality area); AA2: accuracy assessment on the predicted newly emerged mortality 
areas. For training period, the accuracy is tested by the 10-fold cross-validation. AUC: 
area under the curve of a receiver operating characteristic plot. OA: overall accuracy. 
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Projections of future MPB mortality made with the GCM projections suggest that the 
MPB outbreak in Grand County would continue to spread until around year 2015 
(Figure 3, Figure SI). The outbreak area predicted under the RCP8.5 scenario was 
consistently higher than that under RCP4.5. From 2015 to 2050, less variation existed in 
the projections among the GCMs and climate projections, as shown from the smaller 
difference between the upper and lower quartiles. 

5 Discussion 

Our Landsat-derived disturbance maps documented an evolving pattern where a severe 
MPB outbreak emerged in the early 2000s, spread throughout the decade with the 
population reaching epidemic levels and peak mortality around 2006. Based on that 
historical record of forest disturbance (Liang et al., 2014a), and the models developed in 
this study, our research provided a means of simulating the landscape-level outbreak 
pattern over time. The major controls on the observed patterns of MPB mortality during 
the time period of this study included residential housing density, density of adjacent 
MPB mortality in previous years, and climate predictors. The diversity of predictor 
variables indicated the complexity in predicting the incidence of MPB mortality, and we 
anticipate that landscape-level outbreak will become more extensive and severe. 

5.1 Detecting spatiotemporal changes in MPB activity and associated uncertainties 

During the course of this study, we experienced challenges in obtaining an ideal dataset 
that accurately depicted the dynamic extent of MPB mortality at a landscape scale 
because of the heterogenous nature of MPB outbreaks. Mortality is rarely complete 
within forest plots, stands, and/or pixel-level satellite observations. Healthy trees are 
often found next to attacked trees and attacked trees can be present in various stages of 
mortality (e.g., red, gray). In addition to those issues, data with high levels of omission 
and/or commission errors should be avoided for use in descriptive and predictive models 
as they may result in false inferences being made about the underlying mechanisms 
influencing patterns of MPB mortality and result in erroneous predictions. However, 
there are very few publicly available data describing the spatial and temporal patterns of 
MPB mortality, and most studies predicting patterns of beetle occurrence collected their 
response variables via field work, ADS, or image interpretation. Among these, ADS has 
been used more widely because of its availability and information richness, such as host 
species and type of disturbance agent (Hicke et al., 2013, Meddens et al., 2012, Preisler 
et al., 2012; Strohm et al., 2013), but the subjective nature and the limited spatiotemporal 
extent of the ADS data made their integration in our analysis problematic (Hicke et al., 
2012a). Johnson & Ross (2008) suggested that the accuracy of ADS data is most 
acceptable for coarse-scale (>500 m) studies, and less suitable at intennediate scales 
(>50 m), and should be cautiously used at fine spatial scales. We also lacked enough 
time-series samples in the disturbed areas to conduct a quantitative evaluation of the ADS 
data. Nonetheless, our visual examination suggested that there are substantial 
duplications in areas of mortality sketch mapped among different years and that the fine- 
scale heterogeneous patterns of MPB outbreaks are not well characterized. Thus, future 
studies that utilize the ADS data should incorporate methods to assess their uncertainties 
and the impacts they have on analyses relying upon them. 
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Instead of relying on the ADS data, we implemented our own algorithm to track and 
identify forest disturbances with Landsat data (Liang et ah, 2014a). Our method built 
upon the LandTrendr algorithm (Kennedy et ah, 2010), but incorporated methods to 
automatically label changed areas identified by LandTrendr with the type of disturbance 
causing the change. The change detection analysis disturbance maps used in this study 
showed promising results, especially in their ability to produce long-term time series of 
insect mortality. Our next steps will pursue testing and application of the change 
detection analysis for disturbance mapping at greater spatial extents. There have been 
several efforts in providing the broad patterns of forest dynamics across the continent, 
e.g., the North American Forest Dynamics project (Masek et ah, 2013), and a recent 30-m 
resolution global forest change product (Hansen et ah, 2013), but the lack of information 
identifying specific disturbance types limits their usability in studies identifying the 
drivers behind the change. Meanwhile, errors from our mapping procedure may have 
affected the accuracy and applicability of the subsequent models. For instance, in the 
accuracy assessment, we found that although the overall accuracy was as high as 90%, 
the commission errors of MPB mortality were higher than omission errors by an average 
of 10%, and higher omission errors were found in the clearcut land cover type (Liang 
etal., 2014a). This might have resulted in sampling errors that were propagated in the 
GLM results. 

5.2 Driving factors of the dynamic beetle infestation pattern 

Residential housing density, beetle pressure and climate were key predictors in the final 
model after stepwise selection, which reflects the effects of human impacts, biological 
dispersal and physical environmental factors in MPB outbreaks. Except for housing 
density, the function of the other key variables in driving beetle outbreaks has already 
been highlighted in some previous studies. For instance, Preisler etal. (2012) found 
beetle pressure, minimum winter temperature, and two-year cumulative precipitation to 
be important predictors of MPB-caused tree mortality. The positive effect of residential 
housing density on MPB mortality indicated that anthropogenic influences provide a 
positive feedback to beetle outbreaks. We suspect that the positive association exists for 
two reasons. First, tree removal resulting from hazard mitigation, timber harvesting, and 
recreational facility construction such as ski resorts is common in Grand County. 
Increasing fragmentation of remaining forests, leaving them with higher edge-to-area 
ratios (Raffa et ah, 2008) and drier conditions because of higher levels of solar radiation 
(Bone et ah, 2013). Consequently, those forests are more likely to be exposed to mass 
attack by MPB. Second, modem urban environments have greatly altered soil physical 
and biochemical properties and heavier pollutant loads, all of which could hamper tree 
growth and make them more susceptible to insect attacks (Bone et ah, 2013). 

The importance of the number of adjacent pixels with MPB mortality in the previous year 
in predicting MPB mortality in our models is in agreement with the current understanding 
of MPB population dynamics (Aukema et ah, 2008; Walter & Platt, 2013). Greater MPB 
densities allow for mass attack and increase the likelihood of tree mortality regardless of 
the vigor and defense system of host trees. Simard etal. (2012) demonstrated that the 
amount of beetle-killed forest in adjacent areas was a key predictor of subsequent 
mortality, and that beetle density is also a potentially factor limiting stand-scale outbreaks 
from developing into landscape-scale outbreaks (Raffa et ah, 2008). Meanwhile, we 
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observed that among all the eight neighborhood mortality variables, only nm_30m was 
retained in the final GLM model, which indicated to us that short-distance dispersal was 
the dominant mode of expansion of the MPB outbreak in Grand County. Long-distance 
dispersal has been suggested to be crucial in the initiation and early stages of infestations, 
but short-distance dispersal dominates the stage when infestations intensify and 
populations reach epidemic levels (Chen & Walton, 2011). Our starting year was 2000, at 
which time MPB had already formed two outbreak clusters in the northern and southern 
comers of our study area. Because of the relatively important role of the short-distance 
dispersal, it is not surprising that we observed and predicted a continuously expanding 
pattern of MPB mortality in subsequent years instead of isolated MPB infestations. 

Our study also identified a negative relationship with summer precipitation and a positive 
relationship with winter temperature and MPB mortality. These relationships have also 
been found by other studies; higher temperatures foster outbreaks whereas lower 
temperatures depress beetle populations (Bentz et al., 2010; Kaufmann et al., 2008; Raffa 
et al., 2008). A warmer climate will reduce the cold-induced mortality in the adult and 
larval stages, and will accelerate the developmental timing within one generation. The 
thennal changes can also detennine the abundance of the fungal species vectored by 
MPB (Six & Bentz, 2007), which will ultimately affect the success of MPB populations. 
The association between precipitation and MPB mortality is less understood than the 
relationships with temperature. Preisler et al. (2012) found that precipitation in the 
previous year increased the odds of outbreak intensification, which could be related to the 
increased beetle brood production within thicker phloem. In contrast, longer drought 
periods may lead to increased host susceptibility and thus result in a higher probability of 
outbreak intensification. The negative effect of summer precipitation as indicated from 
our study supports the latter statement. 

6 Conclusions 

Evidence has been accumulating to document the contribution of climate change to recent 
increases in the frequency, duration, extent, and severity of insect disturbances (Kurz 
et al., 2008). Despite uncertainty in downscaled forecasts of future climate elements, both 
the RCP 4.5 and 8.5 climate-change scenarios indicated an expanding extent of MPB 
mortality, with the implication that future climate conditions in Grand County, Colorado 
will be more suitable for MPB survival. Little variation existed in the projected area of 
MPB mortality in Grand County between years 2015 and 2050, because a limited amount 
of healthy forest with MPB host species may remain at that point. Because of this, our 
projections for Grand County do not fully depict the situation in the Southern Rocky 
Mountains ecoregion, as there is a substantial amount of lodgepole and other species of 
pine forests that currently remain healthy outside of Grand County. Given that projected 
climate conditions within the ecoregion are likely to follow those within Grand County, 
we anticipate the area of forest with MPB mortality within the ecoregion will increase. 
The influence of climate and weather factors on the beetle-caused tree mortality varies 
among different locations (Creeden et al., 2014), and not all areas across western North 
America are expected to be more suitable to MPB survival as temperatures increase 
(Hicke et al., 2006). Thus, the models developed, and the conclusions drawn from this 
study might not be applicable to other areas. However, our overall approach is applicable 
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to other regions experiencing similar insect outbreaks and can aid in generating consistent 
and high-temporal frequency data on insect mortality and other disturbances impacting 
carbon cycling and other ecosystem services. 

Acknowledgments 

This research was supported by the U.S. Geological Survey, Climate and Land Use 
Mission Area Land Change Science Program (grant number G12AC20085), and 
a National High Technology Grant from China (2009AA12200101). We also would like 
to acknowledge the World Climate Research Programme's Working Group on Coupled 
Modelling, which is responsible for CMIP, and we thank the climate modeling groups 
(listed in Table 2 of this paper) for producing and making available their model output. 
For CMIP, the U.S. Department of Energy's Program for Climate Model Diagnosis and 
Intercomparison provides coordinating support and led development of software 
infrastructure in partnership with the Global Organization for Earth System Science 
Portals. Two anonymous reviewers, Diane Stephens and Julie Ann Beston provided 
insightful comments on a previous draft of this manuscript and their comments helped 
greatly to improve the completeness and clarity of the manuscript. Any use of trade, firm, 
or product names is for descriptive purposes only and does not imply endorsement by the 
U.S. Government. 


80 



2012 


2013 


2014 


2015 


2050 



BCC-CSMI. 


CanESM2 


81 













































Figure SI. The spatial explicit maps of predicted MPB mortality under the RCP8.5 
scenario using the final general linear model and the fourteen global climate model 
outputs. 
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Chapter 5 Summary and broader impacts 

In this dissertation, I incorporated multi-disciplinary knowledge and techniques, 
including remote sensing, geographic information science, forestry, ecology, geostatistics, 
statistics and climate science, to answer some critical questions with regard to the 
mountain pine beetle disturbance in the highly landscape heterogeneous Southern Rocky 
Mountain Ecoregion: 

1) How can satellite images with moderate resolution be applied in monitoring 
the chronic mountain pine beetle disturbance over a decadal long period with 
high accuracy and efficiency? 

2) What biotic and abiotic drivers promote the extensive development and 
progressive MPB outbreaks? And what is the future habitat of mountain pine 
beetle under climate change scenarios? 

3) What are the characteristics of disturbance regime, their interactions and the 
successional pathway in this region? 

To address the first question, I investigated the benefits of trajectory-based classification 
for time-series forest disturbance mapping at the Grand County, Colorado, where a high 
concentration of lodgepole pines and MPB mortality are located. Unlike using pixels as 
the basic unit in the traditional classifications, this method employed the temporal 
segment component derived from the image spectral trajectory. Several advantages were 
presented in this method. Firstly, it is suitable for applications of which images are not 
well radiometrically corrected, since it is designed to remove spikes and enhance the 
signal-to-noise ratio of the often noisy time series. Secondly, because one segment 
represents a same type of event happening over the duration, the illogical land cover 
changes can be minimized. Thirdly, a smaller amount of training samples was needed in 
order to achieve the same accuracy as the single-date classification. We further advanced 
this trajectory based procedure to ecoregion scale, and more efforts were put on 
improving the algorithm efficiency by replacing decision tree labeling procedure with an 
automatic random forest step. A set of segment features and their spatial-temporal 
neighborhood features were fed to the random forest classifier to enhance its 
discrimination ability of abrupt and chronic disturbance types. Another remarkable 
improvement is the application of an effective scene overlapping area sampling strategy, 
by which substantial time and labor efforts can be saved in selecting training data. 

Based on the satellite measures of disturbance and recovery pattern, I was able to 
quantify the relative influence of various anthropogenic, biologic and physical drivers 
that have shaped the spatio-temporal patterns of the MPB outbreaks via an ecological 
niche modeling framework. After the removal of spatial autocorrelation effects and non¬ 
significant predictors, the final model suggested that beetle neighborhood mortality, 
winter mean temperature anomaly, and residential housing density were positively 
associated with MPB mortality, whereas summer precipitation was negatively correlated. 
The four variables reflected the effects of human impacts, biological dispersal and 
physical environmental factors separately in outbreaks. I further identified the forested 
area with high risk for future MPB mortality in Grand County to provide informative 
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future distribution information for management purposes. It was found that the extent of 
MPB would likely to expand under both two future climate change projections. Those 
findings imply that future climate conditions will be more suitable for MPB survival. 

The highly spatially and temporal variable disturbance and recovery pattern as revealed 
from the remote sensing mapping results highlighted the importance of resolving the 
complex relationship between different disturbance events, and their compounded effects 
on the resilience mechanism of pine forest ecosystem. Besides bark beetle, fire is another 
major forest disturbance agent in this region whose occurrence is largely determined by 
the environmental conditions. Although it is presumed that beetle infestation could 
increase the risks of fire hazard, I found that fire severity was largely unrelated to pre-fire 
MPB outbreak severity at the landscape level. Beetles can affect the forest ecosystem in a 
variety of ways, and thus it is hard to predict the synergic effects on fire potential since 
the positive and negative impacts could counteract with each other. Fire severity, 
nevertheless, was found to be negatively correlated with the post-fire beetle severity, and 
can result in milder outbreaks. I further assessed the impacts of singular and multiple 
disturbances on the forest regeneration pattern. Stand-replacing events, including fire and 
clearcut, have a high seedling reestablishment rate, whereas only a small portion of the 
MPB disturbed forests was on a trajectory to full recovery. Interactions between 
processes can amplify or mute the overall effects of changes. In this ecoregion, beetle can 
synergistic with fire in increasing the post-fire regeneration rate of serotinous lodgepole 
pines, which suggests that the beetle disturbance created novel characteristics on the 
forest resilience response to fire. 

This dissertation can be offered as a persuasive and informative case study that 
demonstrates and further expands the power of remote sensing application in monitoring 
and modeling landscape ecosystem changes. The outcomes of this dissertation, from the 
aspect of methodology, map product or ecological findings, could be applied to a broad 
range of scientific domains, such as biological conservation, carbon science, hydrology, 
soil science, and ecology. Below I listed two areas that can directly benefit from this 
research. Besides the scientific merits, the results can also provide meaningful and timely 
information to decision makers, which is especially beneficial for vulnerable zones that 
are facing dramatic environment transformations or climate changes. 

Dynamic land cover mapping and phonological metrics. The automatic trajectory-based 
time-series classification working flow proposed in this dissertation could contribute to 
future fine resolution dynamic land cover mapping, which will involve various spatial 
scales and more diverse thematic classes, and thus higher complexity in classification. 
Dynamic land cover is essential in assessing whether a fluctuation caused by human 
intervention or climate change in the system is within the range of natural variation or not, 
and is also a critical land surface parameter in global earth system models. Numerous 
techniques of land cover change detection have been developed, but no single technique 
works well in all contexts, so does this algorithm as its limitations had been summarized 
previously. Besides conducting time-series mapping as a routine task, our trajectory- 
based classification framework is also capable of providing phonological metrics of all 
ecosystems exhibiting identifiable seasonal or annual phenology, such as the onset, 
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duration and magnitude of greenness increase and decrease, and biomass maximum and 
minimum. Those metrics are valuable in enhancing understanding the landscape 
dynamics over time and improving models of phonological responses to climate change. 

More accurate estimation on carbon and biomass. Accurate depiction of forest 
disturbance and recovery over large extent and continuous time is essential in monitoring 
forest carbon fluxes and understanding the role of forest ecosystems in the global climate. 
Since terrestrial carbon fluxes cannot realistically be measured at a regional or global 
scale using modem techniques such as flux tower, remote sensing offers the most feasible 
way in providing such information. There have been successful cases that integrated 
maps of harvest disturbance history and age derived from Landsat time-series with 
ecosystem models (eg. Cohen et ah, 1996; Masek & Collatz, 2006). However, unlike 
physical disturbances such as fire and windthrow, biotic agent caused mortality is never 
instantaneous under any pathogen and insect pathway, and this inherent feature makes the 
long-term growth reduction due to forest insects absent from ecosystem process models, 
which further lead to the systematic overestimation of carbon uptake and storage across 
wide regions. The present MPB mortality maps generated from this dissertation that are 
superior in the fine spatial resolution (30 m) and high temporal frequency (1 year), could 
be utilized as inputs to the processed-based forest productivity models or terrestrial 
ecosystem models, to provide greater realism to modeled impacts. Additionally, the 
potential range shifts of MPB associated with changing climate as estimated by the 
bioclimatic niche model I developed in Chapter four can also be fed into the ecosystem 
models to explore future trends of the impacts on carbon. 
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